Introduction

# rmarkdown::render("MSH3aso.Rmd",
#                   output_file = file.path(out_dir, "MSH3aso.html"))

Setup environment

# Set the working directory where all WGCNA functions are stored
setwd(dirname(rstudioapi::getActiveDocumentContext()$path))
getwd()
## [1] "/Users/michaelflower/Library/Mobile Documents/com~apple~CloudDocs/Documents/ACL/Collaborations/Emma Bunting/2024.07.03 MSH3 ASO paper plots"

Load packages.

# Required packages
packages <- c("plyr", "dplyr", "ggplot2", "data.table", "readxl", "tidyr", "tidyverse", "janitor", "patchwork",
              "ggdark", "xlsx", "openxlsx", "ggpubr", "writexl", "stringr", "ggpmisc", "drc", "mgcv", "scales",
              "rstatix", "kableExtra", "DT", "zoo", "emmeans", "lme4")

# # Install R packages if required
# install.packages(setdiff(packages, rownames(installed.packages())))

# # Install bioconductor packages if required
# BiocManager::install(setdiff(packages, rownames(installed.packages())))
# invisible(lapply(packages, function(x) library(x, character.only=TRUE)))

# Load
lapply(packages, library, character.only = TRUE)
rm(packages)

## Clear objects
# rm(list = ls()) # This would cause rmarkdown to error
# rm(list = ls()[!grepl("_env$", ls())]) # Remove all objects except environment imports

Set variables

# Data path
data_path_spreadsheet <- "./data/MSH3aso_data.xlsx"
data_path_juliet <- "/Users/michaelflower/Library/Mobile Documents/com~apple~CloudDocs/Documents/ACL/Collaborations/Emma Bunting/2024.05.01 Emma MSH3 ASO FA reanalysis/results/02-juliet.RData"


# Output directory
out_dir <- "./results"


# Plot settings
element_text_size = 22
point_size = 3
hide_ns = TRUE


# Genes
gene_colour <- c("MSH2" = "red", # "#F8766D"
                 "MSH3" = "blue", # "#619CFF"
                 "MSH6" = "#00BA38")
gene_order <- c("MSH2", "MSH3", "MSH6")


# Treatments
treatment_rename <- c("Baseline" = "Baseline",
                      "Vehicle" = "Vehicle",
                      "MSH3aso-0.022uM" = "MSH3 ASO 0.022 µM",
                      "MSH3aso-0.26uM" = "MSH3 ASO 0.26 µM",
                      "MSH3aso-3uM" = "MSH3 ASO 3 µM",
                      "SCRaso-3uM" = "SCR ASO 3 µM",
                      "H2O2-500uM" = "H2O2 500 µM",
                      "SCRaso" = "SCR ASO",
                      "MSH3aso_1161149" = "MSH3 ASO 1161149",
                      "MSH3aso_1161173" = "MSH3 ASO 1161173",
                      "MSH3aso_1161329" = "MSH3 ASO 1161329")
                      
treatment_order <- c("Baseline", "Vehicle", "MSH3aso-0.022uM", "MSH3aso-0.26uM", 
                     "MSH3aso-3uM", "SCRaso-3uM", "SCRaso", "H2O2-500uM",
                     "MSH3aso_1161149", "MSH3aso_1161173", "MSH3aso_1161329")

treatment_colour <- c("Baseline" = "darkgrey",
                      "Vehicle" = "blue",
                      "MSH3aso-0.022uM" = "red",
                      "MSH3aso-0.26uM" = "purple",
                      "MSH3aso-3uM" = "orange",
                      "SCRaso-3uM" = "darkgreen",
                      "H2O2-500uM" = "cyan",
                      "SCRaso" = "darkgreen",
                      "MSH3aso_1161149" = "purple",
                      "MSH3aso_1161173" = "purple",
                      "MSH3aso_1161329" = "purple")

treatment_rename_subscript <- c("Baseline" = "Baseline",
                                "Vehicle" = "Vehicle",
                                "MSH3aso-0.022uM" = "MSH3~ASO~0.022~µM",
                                "MSH3aso-0.26uM" = "MSH3~ASO~0.26~µM",
                                "MSH3aso-3uM" = "MSH3~ASO~3~µM",
                                "SCRaso-3uM" = "SCR~ASO~3~µM",
                                "H2O2-500uM" = "H[2]*O[2]~500~µM")


# Genotypes
genotype_rename <- c("Unedited" = "Unedited",
                     "MSH3ko" = "MSH3-/-",
                     "FAN1ko" = "FAN1-/-")
genotype_order <- c("Unedited", "MSH3ko", "FAN1ko")
genotype_colour <- c("Unedited" = "blue",
                     "MSH3ko" = "darkred",
                     "FAN1ko" = "green")


# Ionis dose
treatment_ionis_rename <- c("Vehicle" = "Vehicle",
                            "3ug" = "MSH3 ASO 3 µg",
                            "10ug" = "MSH3 ASO 10 µg",
                            "30ug" = "MSH3 ASO 30 µg",
                            "SCRaso" = "SCR ASO")
treatment_ionis_order <- c("Vehicle", "3ug", "10ug", "30ug")
treatment_ionis_colour <- c("Vehicle" = "blue",
                            "3ug" = "red",
                            "10ug" = "purple",
                            "30ug" = "orange",
                            "SCRaso" = "darkgreen")


# Tissue
tissue_order <- c("Cortex", "Striatum", "Brainstem", "Spinal cord")

# Output formats
output_formats <- "svg" # c("png", "svg", "eps")

# Models to fit
my_models <- list(
  "linear" = y ~ x,
  "log" = y ~ log(x), # Can't use because instability rate goes negative at low MSH3 levels
  "poly2" = y ~ poly(x, 2),
  # "poly2raw" = y ~ poly(x, 2, raw = TRUE),
  "poly3" = y ~ poly(x, 3),
  # "poly3raw" = y ~ poly(x, 3, raw = TRUE),
  "exponential" = I(log(y)) ~ x,  # Exponential model # Can't use because instability rate goes negative at low MSH3 levels
  # "logistic" = I(log(y / (1 - y))) ~ x,  # Simple logistic model # Can't use because instability rate goes negative at low MSH3 levels
  "exponential" = y ~ exp(x)
)

# Log increment
log_increment = 0.001

# p adjust method
padj_method = "BH" # "bonferroni"

# Errorbar metric
errorbar_metric = "SE" # CL, SD

Import data

Spreadsheet

# Function to import all sheets
read_full_excel <- function(file_path = NULL){ 
  
  # Imports the name of the sheets in the excel file and creates a vector
  sheet_names <- readxl::excel_sheets(file_path)
  
  # Create list and add each excel sheet as a dataframe by a loop
  sheets <- list()
  for(i in 1:length(sheet_names)){
    sheets[[i]] <-  readxl::read_xlsx(file_path, sheet = sheet_names[i])
  }
  
  # Add the sheet names to the list: 
  names(sheets) <- sheet_names
  
  return(sheets)
  
}

# Import
MSH3aso_data <- read_full_excel(file_path = data_path_spreadsheet)

Juliet

if (file.exists(file.path(out_dir, "juliet_minimal.RData"))) {
  load(file.path(out_dir, "juliet_minimal.RData"))
} else {
  
  # Import juliet data
  if(!exists("juliet_env")) {
    juliet_env <- new.env()
    load(file = data_path_juliet,
         envir = juliet_env)
    ls(envir = juliet_env)
  }
  
  # Extract required objects
  timecourse_stats <- juliet_env$timecourse_stats
  timecourse_stats_manual <- juliet_env$timecourse_stats_manual
  juliet_annotation <- juliet_env$annotation
  
  # Extract slopes
  my_stats <- timecourse_stats$predvar2$ii
  my_experment_names <- setNames(names(my_stats), names(my_stats))
  juliet_slopes <- lapply(my_experment_names, function(my_experiment_name,
                                               my_stats) {
    my_stats[[my_experiment_name]]$slopes_uncorrected %>%
      dplyr::mutate(experiment_name = my_experiment_name) %>%
      relocate(experiment_name)
  },
  my_stats = my_stats)
  juliet_slopes <- rbindlist(juliet_slopes)
  
  
  # Extract slopes from manually curated MSH3ko vs CRISPRwt experiment
  my_slopes <- timecourse_stats_manual$predvar2$ii$slopes_uncorrected %>%
    dplyr::mutate(experiment_name = "MSH3ko_CRISPRwt_manual") %>%
    relocate(experiment_name)
  juliet_slopes <- rbind(juliet_slopes, my_slopes)
  
  # Save
  save(juliet_slopes, juliet_annotation,
       file = file.path(out_dir, "juliet_minimal.RData"))
  
}

Model fitting function

Define a function to fit models.

# Vector of model names
model_names <- setNames(names(my_models), names(my_models))

# Function to plot fits
plot_model_function <- function(model_name,
                                models,
                                data,
                                xvar, xlabel,
                                yvar, ylabel,
                                group_var = NA, group_label = NA,
                                log_increment = 0.001,
                                plot_se = TRUE) {

  # # Manual
  # model_name = model_names[[2]]
  # models = my_models
  # data = plot_data
  # xvar = "dose_uM"
  # xlabel = "MSH3 ASO dose (µM)"
  # yvar = "level"
  # ylabel = "Relative protein level (%)"
  # group_var = NA
  # group_label = NA
  # log_increment = 0.001
  # plot_se = TRUE

  # # Report
  # print(model_name)

  # Extract formula
  my_model <- models[[model_name]]

  # If formula contains log, use the pre-logged data with +0.001
  my_model_str <- deparse(my_model)
  contains_log <- grepl("log", my_model_str)
  if (contains_log) {
    data <- data %>%
      dplyr::mutate(!!sym(xvar) := !!sym(xvar) + log_increment)
  }

  # Plot
  if (is.na(group_var)) {

    my_plot <-
      ggplot(data, aes(x = !!sym(xvar), y = !!sym(yvar))) +
      # geom_point(size = point_size) +
      geom_point() +
      stat_smooth(method = lm,
                  formula = my_model,
                  fullrange = TRUE, se = plot_se) +
      stat_poly_eq(formula = my_model,
                   aes(label = paste(after_stat(eq.label),
                                     after_stat(rr.label),
                                     after_stat(p.value.label),
                                     sep = "~~~")),
                   parse = TRUE, coef.digits = 3, f.digits = 3, p.digits = 3,
                   rr.digits = 3, size = 3) +
      labs(title = model_name,
           subtitle = c(my_model),
           x = xlabel,
           y = ylabel) +
      theme_minimal()

  } else {

     my_plot <-
       ggplot(data, aes(x = !!sym(xvar), y = !!sym(yvar), colour = !!sym(group_var), fill = !!sym(group_var))) +
       geom_point() +
       stat_smooth(method = lm,
                   formula = my_model,
                   fullrange = TRUE, se = plot_se) +
       stat_poly_eq(formula = my_model,
                    aes(label = paste(after_stat(eq.label),
                                      after_stat(rr.label),
                                      after_stat(p.value.label),
                                      sep = "~~~")),
                    parse = TRUE, coef.digits = 3, f.digits = 3, p.digits = 3,
                    rr.digits = 3, size = 3) +
       labs(title = model_name,
            subtitle = c(my_model),
            x = xlabel,
            y = ylabel,
            colour = group_label,
            fill = group_label) +
       theme_minimal()

  }

  # Convert xy notation to real variables for lm function
  left <- gsub("\\by\\b", yvar, as.character(my_model[2]))
  right <- gsub("\\bx\\b", xvar, as.character(my_model[3])) # This pattern \\bx\\b ensures that "x" is replaced only when it stands alone as a word. \\b denotes a word boundary.

  # Linear model
  if (is.na(group_var)) {
    my_model_long <- sprintf("%s ~ %s", left, right) # https://stackoverflow.com/questions/26381410/edit-and-reuse-the-formula-part-of-the-call-for-a-model-in-r
    my_model_long <- as.formula(my_model_long)
    my_lm <- lm(my_model_long, data = data)
    r2 <- summary(my_lm)$r.squared
  } else {
    my_model_long <- sprintf("%s ~ %s * %s", left, right, group_var) # https://stackoverflow.com/questions/26381410/edit-and-reuse-the-formula-part-of-the-call-for-a-model-in-r
    my_model_long <- as.formula(my_model_long)
    my_lm <- lm(my_model_long, data = data)
    r2 <- summary(my_lm)$r.squared
  }

  # Output
  out <- mget(c("my_plot", "my_model", "r2"))
  return(out)

}

Titration western MSH3 1wk

Fit models

# Experiment ID
experiment_id <- names(MSH3aso_data)[[1]]
experiment_id
## [1] "titration_western_MSH3_1wk"
# Make results directory
dir.create(file.path(out_dir, experiment_id),
           recursive = TRUE)

# Format data
plot_data <- MSH3aso_data[[experiment_id]] %>%
  pivot_longer(!dose_uM,
               names_to = "label",
               values_to = "level") %>%
  separate_wider_delim(label, "_", names = c("gene", "rep")) %>%
  dplyr::mutate(dose_uM = as.numeric(dose_uM),
                level = as.numeric(level),
                dose_uM_forlog = dose_uM + log_increment) %>% # Add increment (+0.001) to log the data
  dplyr::filter(!is.na(level))

# Plot using each formula
my_fits <- lapply(model_names, plot_model_function,
                  models = my_models,
                  data = plot_data,
                  xvar = "dose_uM",
                  xlabel = "MSH3 ASO dose (µM)",
                  yvar = "level",
                  ylabel = "Relative protein level (%)",
                  group_var = NA,
                  group_label = NA,
                  log_increment = log_increment,
                  plot_se = TRUE)

# Plot model fits
my_plots <- lapply(my_fits, function(my_fit) {
  my_fit$my_plot
})
wrap_plots(my_plots) +
  plot_annotation(title = "Regression models",
                  theme = theme(plot.title = element_text(hjust = 0.5, face = "bold")))

# r2 table
my_r2 <- lapply(my_fits, function(my_fit) {
  my_fit[["r2"]]
})
r2_table <- data.frame(model_name = names(my_r2),
                       r2 = unlist(my_r2)) %>%
  tibble::rownames_to_column("model_number") %>%
  arrange(-r2)

# Display r2 table
kbl(r2_table, caption = "R-squared") %>%
  kable_styling(bootstrap_options = "striped")
R-squared
model_number model_name r2
2 log 0.8776671
4 poly3 0.8705157
5 exponential 0.7160194
6 exponential 0.7160194
3 poly2 0.6906890
1 linear 0.5987664

Functions

# Helper function to calculate the effective dose for a given response level using interpolation
ed_interpolation_helper <- function(data, response_var, pred_var, formula, model_function, response_level) {
  
  if (model_function == "lm") {
    my_model <- lm(formula, data = data)
  } else if (model_function == "drm") {
    my_model <- drm(formula, data = data, fct = LL.4())
  }
  predicted_values <- setNames(data.frame(seq(0, max(data[[pred_var]]), length.out = 1000)), pred_var)
  predictions <- predict(my_model, newdata = predicted_values, interval = "confidence")
  predicted_values <- cbind(predicted_values, predictions)
  colnames(predicted_values) <- c(pred_var, response_var, "lower_ci", "upper_ci")
  ed <- approx(predicted_values[, response_var], predicted_values[[pred_var]], xout = response_level)$y
  return(list(predicted_values = predicted_values,
              ed = ed))
}


# Effective dose bootstrap function to calculate confidence intervals
ed_bootstrap_helper <- function(data, response_var, pred_var, formula, model_function, response_level, n_bootstrap = 1000, confidence_level = 0.95) {
  
  # Resample data and compute ED for each resample
  ed_bs <- replicate(n_bootstrap, {
    data_resample <- data[sample(nrow(data), replace = TRUE), ]
    my_ed <- ed_interpolation_helper(data = data_resample,
                                     response_var = response_var,
                                     pred_var = pred_var,
                                     formula = formula,
                                     model_function = model_function,
                                     response_level = response_level)
    return(my_ed$ed)
  })
  
  # Calculate confidence intervals
  alpha <- 1 - confidence_level
  ci_lower = quantile(na.omit(ed_bs), probs = alpha / 2)
  ci_upper = quantile(na.omit(ed_bs), probs = 1 - (alpha / 2))
  
  # Output
  return(list(ci_lower = as.numeric(ci_lower),
              ci_upper = as.numeric(ci_upper)))
}


# Bring the effective dose (ED) interpolation and confidence interval functions together
ed_interpolation <- function(data, response_var, pred_var, formula, model_function, response_level, n_bootstrap = 1000, confidence_level = 0.95) {
  
  # Calculate effective dose
  ed_result <- ed_interpolation_helper(data = data,
                                       response_var = response_var,
                                       pred_var = pred_var,
                                       formula = formula,
                                       model_function = model_function,
                                       response_level = response_level)
  
  # Confidence intervals
  ed_ci <- ed_bootstrap_helper(data = data,
                               response_var = response_var,
                               pred_var = pred_var,
                               formula = formula,
                               model_function = model_function,
                               response_level = response_level,
                               n_bootstrap = n_bootstrap,
                               confidence_level = confidence_level)
  
  # Output
  return(list(ed = ed_result$ed,
                   ci_lower = ed_ci$ci_lower,
                   ci_upper = ed_ci$ci_upper,
                   predicted_values = ed_result$predicted_values))
}

IC50

# Dose response model
my_drm <- drm(level ~ dose_uM, data = plot_data, fct = LL.4())

# Extract IC50 value
ic50 <- ED(my_drm, 50)
## 
## Estimated effective doses
## 
##        Estimate Std. Error
## e:1:50   1.7779    10.8055
ic50_estimate <- ic50[[1]]
ic50_se <- ic50[[2]]
ic50_upper <- ic50_estimate + 1.96 * ic50_se
ic50_lower <- ic50_estimate - 1.96 * ic50_se

# Plot drm model
plot(my_drm, main = "Dose-response curve using drc")
abline(h = 50, col = "blue", lty = 2)
abline(v = ic50_estimate, col = "blue", lty = 2)

# Calculate IC50 by interpolation
ic50_interpolation <- ed_interpolation(data = plot_data,
                                       response_var = "level",
                                       pred_var = "dose_uM_forlog",
                                       formula = as.formula("level ~ log(dose_uM_forlog)"),
                                       model_function = "lm",
                                       response_level = 50,
                                       n_bootstrap = 1000,
                                       confidence_level = 0.95)

# Readjust for the log increment
ic50_interpolation$ed <- ic50_interpolation$ed - log_increment
ic50_interpolation$ci_lower <- ic50_interpolation$ci_lower - log_increment
ic50_interpolation$ci_upper <- ic50_interpolation$ci_upper - log_increment

# Plot linear model
plot(level ~ dose_uM_forlog, 
     data = ic50_interpolation$predicted_values,
     main = "Dose-response curve using interpolation")
abline(h = 50, col = "blue", lty = 2)
abline(v = ic50_interpolation$ed, col = "blue", lty = 2)
abline(v = ic50_interpolation$ci_lower, col = "lightblue", lty = 2)
abline(v = ic50_interpolation$ci_upper, col = "lightblue", lty = 2)

# Table of IC50
ic50_table <- data.frame(method = c("drc", "lm_interpolation"),
                         IC50 = c(ic50_estimate, ic50_interpolation$ed),
                         ci_lower = c(ic50_lower, ic50_interpolation$ci_lower),
                         ci_upper = c(ic50_upper, ic50_interpolation$ci_upper))

# Display IC50 table
kbl(ic50_table, caption = "IC50") %>%
  kable_styling(bootstrap_options = "striped")
IC50
method IC50 ci_lower ci_upper
drc 1.7779260 -19.4008695 22.9567214
lm_interpolation 0.2108101 0.1296216 0.3853176

Linear dose axis

# Plot
my_plot <-
  ggplot(plot_data, aes(x = dose_uM_forlog, y = level)) +
  geom_point(aes(colour = gene), size = point_size) +
  stat_smooth(aes(colour = gene, fill = gene),
              method = lm, formula = y ~ log(x), fullrange = TRUE, 
              size = 1, alpha = 0.1) +
  geom_segment(aes(x = -Inf, xend = ic50_interpolation$ed, 
                   y = 50, yend = 50), 
               colour = "black", linetype = "dashed", size = 0.5) +
  geom_segment(aes(x = ic50_interpolation$ed, xend = ic50_interpolation$ed, 
                   y = -Inf, yend = 50), 
               colour = "black", linetype = "dashed", size = 0.5) +
  geom_rect(aes(xmin = ic50_interpolation$ci_lower, xmax = ic50_interpolation$ci_upper,
                ymin = -Inf, ymax = 50), 
            fill = "darkgrey", alpha = 0.03, colour = NA) +
  scale_colour_manual(values = gene_colour) +
  scale_fill_manual(values = gene_colour) +
  scale_x_continuous(breaks = pretty_breaks()) +
  scale_y_continuous(breaks = pretty_breaks()) +
  labs(x = "MSH3 ASO dose (µM)",
       y = "Relative MSH3 protein level (%)",
       colour = "Protein",
       fill = "Protein") +
  theme_bw() +
  theme(text = element_text(size = element_text_size))
my_plot

# Export
for (output_format in output_formats) {
  ggsave(filename = file.path(out_dir, experiment_id, paste0(experiment_id, "_linearx.", output_format)),
         plot = my_plot)
}
## Saving 7 x 5 in image
# Plot without legend
my_plot <- my_plot +
  theme(legend.position = "none")
my_plot

# Export
for (output_format in output_formats) {
  ggsave(filename = file.path(out_dir, experiment_id, paste0(experiment_id, "_linearx_nolegend.", output_format)),
         plot = my_plot,
         height = 6, width = 6)
}

Log dose axis

# Plot
my_plot <-
  ggplot(plot_data, aes(x = dose_uM_forlog, y = level, colour = gene, fill = gene)) +
  geom_point(size = point_size) +
  stat_smooth(method = lm, 
              formula = y ~ x, fullrange = TRUE, 
              size = 1, alpha = 0.1) +
  geom_segment(aes(x = 0.001, xend = ic50_interpolation$ed, 
                   y = 50, yend = 50), 
               colour = "black", linetype = "dashed", size = 0.5) +
  geom_segment(aes(x = ic50_interpolation$ed, xend = ic50_interpolation$ed, 
                   y = -Inf, yend = 50), 
               colour = "black", linetype = "dashed", size = 0.5) +
  geom_rect(aes(xmin = ic50_interpolation$ci_lower, xmax = ic50_interpolation$ci_upper,
                ymin = -Inf, ymax = 50), 
            fill = "darkgrey", alpha = 0.03, colour = NA) +
  scale_colour_manual(values = gene_colour) +
  scale_fill_manual(values = gene_colour) +
  scale_x_continuous(trans = "log2",
                     breaks = c(0.001, 0.01, 0.1, 1, 2, 3),
                     labels = function(x) {
                       ifelse(x < 1, round(x, digits = 3), round(x, digits = 0))
                     }) +
  scale_y_continuous(breaks = pretty_breaks()) +
  labs(x = "MSH3 ASO dose (µM)",
       y = "Relative MSH3 protein level (%)",
       colour = "Protein",
       fill = "Protein") +
  theme_bw() +
  theme(text = element_text(size = element_text_size))
my_plot

# Export
for (output_format in output_formats) {
  ggsave(filename = file.path(out_dir, experiment_id, paste0(experiment_id, "_logx.", output_format)),
         plot = my_plot)
}
## Saving 7 x 5 in image
# Plot without legend
my_plot <- my_plot +
  theme(legend.position = "none")
my_plot

# Export
for (output_format in output_formats) {
  ggsave(filename = file.path(out_dir, experiment_id, paste0(experiment_id, "_logx_nolegend.", output_format)),
         plot = my_plot,
         height = 6, width = 6)
}

Titration western MSH2 and MSH6 1wk

Fit models

# Experiment ID
experiment_id <- names(MSH3aso_data)[[2]]
experiment_id
## [1] "titration_western_MSH2_MSH6_1wk"
# Make results directory
dir.create(file.path(out_dir, experiment_id),
           recursive = TRUE)

# Format data
plot_data <- MSH3aso_data[[experiment_id]] %>%
  pivot_longer(!dose_uM,
               names_to = "label",
               values_to = "level") %>%
  separate_wider_delim(label, "_", names = c("gene", "rep")) %>%
  dplyr::mutate(dose_uM = as.numeric(dose_uM),
                level = as.numeric(level),
                dose_uM_forlog = dose_uM + log_increment) %>% # Add increment (+0.001) to log the data
  dplyr::filter(!is.na(level))

# Relevel variables
plot_data <- plot_data %>%
  dplyr::mutate(gene = fct_relevel(gene, intersect(gene_order, unique(plot_data$gene))))

# Plot using each formula
my_fits <- lapply(model_names, plot_model_function,
                  models = my_models,
                  data = plot_data,
                  xvar = "dose_uM",
                  xlabel = "MSH3 ASO dose (µM)",
                  yvar = "level",
                  ylabel = "Relative protein level (%)",
                  group_var = "gene",
                  group_label = "Gene",
                  log_increment = log_increment,
                  plot_se = TRUE)

# Plot model fits
my_plots <- lapply(my_fits, function(my_fit) {
  my_fit$my_plot
})
wrap_plots(my_plots) +
  plot_annotation(title = "Regression models",
                  theme = theme(plot.title = element_text(hjust = 0.5, face = "bold")))

# r2 table
my_r2 <- lapply(my_fits, function(my_fit) {
  my_fit[["r2"]]
})
r2_table <- data.frame(model_name = names(my_r2),
                       r2 = unlist(my_r2)) %>%
  tibble::rownames_to_column("model_number") %>%
  arrange(-r2)

# Display r2 table
kbl(r2_table, caption = "R-squared") %>%
  kable_styling(bootstrap_options = "striped")
R-squared
model_number model_name r2
4 poly3 0.1953677
3 poly2 0.1478138
2 log 0.1371533
5 exponential 0.1025985
6 exponential 0.1025985
1 linear 0.1009314

Linear dose axis

# Plot
my_plot <-
  ggplot(plot_data, aes(x = dose_uM, y = level, colour = gene, fill = gene)) +
  geom_point(size = point_size) +
  stat_smooth(method = lm, 
              formula = y ~ x,
              fullrange = TRUE, alpha = 0.2) +
  scale_colour_manual(values = gene_colour) +
  scale_fill_manual(values = gene_colour) +
  scale_x_continuous(breaks = pretty_breaks()) +
  scale_y_continuous(breaks = pretty_breaks()) +
  labs(x = "MSH3 ASO dose (µM)",
       y = "Relative protein level (%)",
       colour = "Protein",
       fill = "Protein") +
  theme_bw() +
  theme(text = element_text(size = element_text_size))
my_plot

# Export
for (output_format in output_formats) {
  ggsave(filename = file.path(out_dir, experiment_id, paste0(experiment_id, "_linearx.", output_format)),
         plot = my_plot)
}
## Saving 7 x 5 in image
# Plot without legend
my_plot <- my_plot +
  theme(legend.position = "none")
my_plot

# Export
for (output_format in output_formats) {
  ggsave(filename = file.path(out_dir, experiment_id, paste0(experiment_id, "_linearx_nolegend.", output_format)),
         plot = my_plot,
         height = 6, width = 6)
}

Log dose axis

# Plot
my_plot <-
  ggplot(plot_data, aes(x = dose_uM, y = level, colour = gene, fill = gene)) +
  geom_point(size = point_size) +
  stat_smooth(method = lm, 
              formula = y ~ x,
              fullrange = TRUE, alpha = 0.2) +
  scale_colour_manual(values = gene_colour) +
  scale_fill_manual(values = gene_colour) +
  scale_x_continuous(trans = "log2",
                     breaks = c(0.001, 0.01, 0.1, 1, 2, 3),
                     labels = function(x) {
                       ifelse(x < 1, round(x, digits = 3), round(x, digits = 0))
                     }) +
  scale_y_continuous(breaks = pretty_breaks()) +
  labs(x = "MSH3 ASO dose (µM)",
       y = "Relative protein level (%)",
       colour = "Protein",
       fill = "Protein") +
  theme_bw() +
  theme(text = element_text(size = element_text_size))
my_plot

# Export
for (output_format in output_formats) {
  ggsave(filename = file.path(out_dir, experiment_id, paste0(experiment_id, "_logx.", output_format)),
         plot = my_plot)
}
## Saving 7 x 5 in image
# Plot without legend
my_plot <- my_plot +
  theme(legend.position = "none")
my_plot

# Export
for (output_format in output_formats) {
  ggsave(filename = file.path(out_dir, experiment_id, paste0(experiment_id, "_logx_nolegend.", output_format)),
         plot = my_plot,
         height = 6, width = 6)
}

Titration western MSH3 9wk

Fit models

# Experiment ID
experiment_id <- names(MSH3aso_data)[[3]]
experiment_id
## [1] "titration_western_MSH3_9wk"
# Make results directory
dir.create(file.path(out_dir, experiment_id),
           recursive = TRUE)

# Format data
plot_data <- MSH3aso_data[[experiment_id]] %>%
  pivot_longer(!dose_uM,
               names_to = "label",
               values_to = "level") %>%
  separate_wider_delim(label, "_", names = c("gene", "rep")) %>%
  dplyr::mutate(dose_uM = as.numeric(dose_uM),
                level = as.numeric(level),
                dose_uM_forlog = dose_uM + log_increment) %>% # Add increment (+0.001) to log the data
  dplyr::filter(!is.na(level))

# Plot using each formula
my_fits <- lapply(model_names, plot_model_function,
                  models = my_models,
                  data = plot_data,
                  xvar = "dose_uM",
                  xlabel = "MSH3 ASO dose (µM)",
                  yvar = "level",
                  ylabel = "Relative protein level (%)",
                  group_var = NA,
                  group_label = NA,
                  log_increment = log_increment,
                  plot_se = TRUE)

# Plot model fits
my_plots <- lapply(my_fits, function(my_fit) {
  my_fit$my_plot
})
wrap_plots(my_plots) +
  plot_annotation(title = "Regression models",
                  theme = theme(plot.title = element_text(hjust = 0.5, face = "bold")))

# r2 table
my_r2 <- lapply(my_fits, function(my_fit) {
  my_fit[["r2"]]
})
r2_table <- data.frame(model_name = names(my_r2),
                       r2 = unlist(my_r2)) %>%
  tibble::rownames_to_column("model_number") %>%
  arrange(-r2)

# Display r2 table
kbl(r2_table, caption = "R-squared") %>%
  kable_styling(bootstrap_options = "striped")
R-squared
model_number model_name r2
4 poly3 0.9217185
2 log 0.9204592
3 poly2 0.8405180
5 exponential 0.7286988
6 exponential 0.7286988
1 linear 0.5735873

IC50

# Dose response model
my_drm <- drm(level ~ dose_uM, data = plot_data, fct = LL.4())

# Extract IC50 value
ic50 <- ED(my_drm, 50)
## 
## Estimated effective doses
## 
##        Estimate Std. Error
## e:1:50  0.13746    0.35855
ic50_estimate <- ic50[[1]]
ic50_se <- ic50[[2]]
ic50_upper <- ic50_estimate + 1.96 * ic50_se
ic50_lower <- ic50_estimate - 1.96 * ic50_se

# Plot drm model
plot(my_drm, main = "Dose-response curve using drc")
abline(h = 50, col = "green", lty = 2)
abline(v = ic50_estimate, col = "blue", lty = 2)

# Calculate IC50 by interpolation
ic50_interpolation <- ed_interpolation(data = plot_data,
                                       response_var = "level",
                                       pred_var = "dose_uM_forlog",
                                       formula = as.formula("level ~ log(dose_uM_forlog)"),
                                       model_function = "lm",
                                       response_level = 50,
                                       n_bootstrap = 1000,
                                       confidence_level = 0.95)

# Readjust for the log increment
ic50_interpolation$ed <- ic50_interpolation$ed - log_increment
ic50_interpolation$ci_lower <- ic50_interpolation$ci_lower - log_increment
ic50_interpolation$ci_upper <- ic50_interpolation$ci_upper - log_increment

# Plot linear model
plot(level ~ dose_uM_forlog, 
     data = ic50_interpolation$predicted_values,
     main = "Dose-response curve using interpolation")
abline(h = 50, col = "blue", lty = 2)
abline(v = ic50_interpolation$ed, col = "blue", lty = 2)
abline(v = ic50_interpolation$ci_lower, col = "lightblue", lty = 2)
abline(v = ic50_interpolation$ci_upper, col = "lightblue", lty = 2)

# Table of IC50
ic50_table <- data.frame(method = c("drc", "lm_interpolation"),
                         IC50 = c(ic50_estimate, ic50_interpolation$ed),
                         ci_lower = c(ic50_lower, ic50_interpolation$ci_lower),
                         ci_upper = c(ic50_upper, ic50_interpolation$ci_upper))

# Display IC50 table
kbl(ic50_table, caption = "IC50") %>%
  kable_styling(bootstrap_options = "striped")
IC50
method IC50 ci_lower ci_upper
drc 0.1374605 -0.5652906 0.8402115
lm_interpolation 0.0915830 0.0436033 0.1893050

Linear dose axis

# Plot
my_plot <-
  ggplot(plot_data, aes(x = dose_uM_forlog, y = level, colour = gene, fill = gene)) +
  geom_point(size = point_size) +
  stat_smooth(method = lm, 
              formula = y ~ log(x), fullrange = TRUE, 
              size = 1, alpha = 0.1) +
  geom_segment(aes(x = -Inf, xend = ic50_interpolation$ed, 
                   y = 50, yend = 50), 
               colour = "black", linetype = "dashed", size = 0.5) +
  geom_segment(aes(x = ic50_interpolation$ed, xend = ic50_interpolation$ed, 
                   y = -Inf, yend = 50), 
               colour = "black", linetype = "dashed", size = 0.5) +
  geom_rect(aes(xmin = ic50_interpolation$ci_lower, xmax = ic50_interpolation$ci_upper,
                ymin = -Inf, ymax = 50), 
            fill = "darkgrey", alpha = 0.03, colour = NA) +
  
  
  # geom_rect(aes(xmin = ic50_interpolation$ci_lower, xmax = ic50_interpolation$ci_upper,
  #               ymin = 0, ymax = Inf), fill = "grey", alpha = 0.05, colour = NA) +
  # geom_hline(yintercept = 50, colour = "black", linetype = "dashed") +
  # geom_vline(xintercept = ic50_interpolation$ed, colour = "black", linetype = "dashed") +
  scale_colour_manual(values = gene_colour) +
  scale_fill_manual(values = gene_colour) +
  scale_x_continuous(breaks = pretty_breaks()) +
  scale_y_continuous(breaks = pretty_breaks()) +
  labs(x = "MSH3 ASO dose (µM)",
       y = "Relative MSH3 protein level (%)",
       colour = "Protein",
       fill = "Protein") +
  theme_bw() +
  theme(text = element_text(size = element_text_size))
my_plot

# Export
for (output_format in output_formats) {
  ggsave(filename = file.path(out_dir, experiment_id, paste0(experiment_id, "_linearx.", output_format)),
         plot = my_plot)
}
## Saving 7 x 5 in image
# Plot without legend
my_plot <- my_plot +
  theme(legend.position = "none")
my_plot

# Export
for (output_format in output_formats) {
  ggsave(filename = file.path(out_dir, experiment_id, paste0(experiment_id, "_linearx_nolegend.", output_format)),
         plot = my_plot,
         height = 6, width = 6)
}

Log dose axis

# Plot
my_plot <-
  ggplot(plot_data, aes(x = dose_uM_forlog, y = level, colour = gene, fill = gene)) +
  geom_point(size = point_size) +
  stat_smooth(method = lm, 
              formula = y ~ x, fullrange = TRUE, 
              size = 1, alpha = 0.1) +
  geom_segment(aes(x = 0.001, xend = ic50_interpolation$ed, 
                   y = 50, yend = 50), 
               colour = "black", linetype = "dashed", size = 0.5) +
  geom_segment(aes(x = ic50_interpolation$ed, xend = ic50_interpolation$ed, 
                   y = -Inf, yend = 50), 
               colour = "black", linetype = "dashed", size = 0.5) +
  geom_rect(aes(xmin = ic50_interpolation$ci_lower, xmax = ic50_interpolation$ci_upper,
                ymin = -Inf, ymax = 50), 
            fill = "darkgrey", alpha = 0.03, colour = NA) +
  scale_colour_manual(values = gene_colour) +
  scale_fill_manual(values = gene_colour) +
  scale_x_continuous(trans = "log2",
                     breaks = c(0.001, 0.01, 0.1, 1, 2, 3),
                     labels = function(x) {
                       ifelse(x < 1, round(x, digits = 3), round(x, digits = 0))
                     }) +
  scale_y_continuous(breaks = pretty_breaks()) +
  labs(x = "MSH3 ASO dose (µM)",
       y = "Relative MSH3 protein level (%)",
       colour = "Protein",
       fill = "Protein") +
  theme_bw() +
  theme(text = element_text(size = element_text_size))
my_plot

#  I'M HERE!

# Export
for (output_format in output_formats) {
  ggsave(filename = file.path(out_dir, experiment_id, paste0(experiment_id, "_logx.", output_format)),
         plot = my_plot)
}
## Saving 7 x 5 in image
# Plot without legend
my_plot <- my_plot +
  theme(legend.position = "none")
my_plot

# Export
for (output_format in output_formats) {
  ggsave(filename = file.path(out_dir, experiment_id, paste0(experiment_id, "_logx_nolegend.", output_format)),
         plot = my_plot,
         height = 6, width = 6)
}

Duration response

Fit models

# Experiment ID
experiment_id <- names(MSH3aso_data)[[4]]
experiment_id
## [1] "duration_response"
# Make results directory
dir.create(file.path(out_dir, experiment_id),
           recursive = TRUE)

# Format data
plot_data <- MSH3aso_data[[experiment_id]] %>%
  pivot_longer(!treatment_duration_weeks,
               names_to = "label",
               values_to = "level") %>%
  separate_wider_delim(label, "_", names = c("treatment", "rep")) %>%
  dplyr::mutate(treatment_duration_weeks = as.numeric(treatment_duration_weeks),
                level = as.numeric(level) * 100) %>%
  dplyr::filter(!is.na(level))

# Relevel variables
plot_data <- plot_data %>%
  dplyr::mutate(treatment = fct_relevel(treatment, intersect(treatment_order, unique(plot_data$treatment))))


# Plot using each formula
my_fits <- lapply(model_names, plot_model_function,
                  models = my_models,
                  data = plot_data,
                  xvar = "treatment_duration_weeks",
                  xlabel = "Week",
                  yvar = "level",
                  ylabel = "Relative expression (%)",
                  group_var = "treatment",
                  group_label = "Treatment",
                  log_increment = log_increment,
                  plot_se = TRUE)

# Plot model fits
my_plots <- lapply(my_fits, function(my_fit) {
  my_fit$my_plot
})
wrap_plots(my_plots) +
  plot_annotation(title = "Regression models",
                  theme = theme(plot.title = element_text(hjust = 0.5, face = "bold")))

# r2 table
my_r2 <- lapply(my_fits, function(my_fit) {
  my_fit[["r2"]]
})
r2_table <- data.frame(model_name = names(my_r2),
                       r2 = unlist(my_r2)) %>%
  tibble::rownames_to_column("model_number") %>%
  arrange(-r2)

# Display r2 table
kbl(r2_table, caption = "R-squared") %>%
  kable_styling(bootstrap_options = "striped")
R-squared
model_number model_name r2
4 poly3 0.7813014
3 poly2 0.7771263
2 log 0.7704015
1 linear 0.7668544
5 exponential 0.7420702
6 exponential 0.7420702

Plot

# Plot
my_plot <-
  ggplot(plot_data, aes(x = treatment_duration_weeks, y = level, colour = treatment, fill = treatment)) +
  geom_point(size = point_size) +
  stat_smooth(method = lm, 
              formula = y ~ log(x),
              fullrange = FALSE, alpha = 0.2) +
  scale_colour_manual(values = treatment_colour,
                      labels = treatment_rename) +
  scale_fill_manual(values = treatment_colour,
                      labels = treatment_rename) +
  scale_x_continuous(breaks = c(0, 3, 6, 9, 12, 15)) +
  scale_y_continuous(breaks = pretty_breaks()) +
  labs(x = "Weeks",
       y = "Relative expression (%)",
       colour = "Treatment",
       fill = "Treatment") +
  theme_bw() +
  theme(text = element_text(size = element_text_size))
my_plot

# Export
for (output_format in output_formats) {
  ggsave(filename = file.path(out_dir, experiment_id, paste0(experiment_id, ".", output_format)),
         plot = my_plot)
}
## Saving 10 x 7 in image
# Plot without legend
my_plot <- my_plot +
  theme(legend.position = "none")
my_plot

# Export
for (output_format in output_formats) {
  ggsave(filename = file.path(out_dir, experiment_id, paste0(experiment_id, "_nolegend.", output_format)),
         plot = my_plot,
         height = 6, width = 6)
}

Statistics

Descriptive statistics

Mean, SD and SEM MSH3 level at each timepoint, for each treatment.

# Descriptive statistics
descriptive_stats <- plot_data %>%
  group_by(treatment_duration_weeks, treatment) %>%
  summarise(mean = mean(level),
            sd = sd(level),
            sem = sd(level)/sqrt(n())) %>%
  ungroup()
## `summarise()` has grouped output by 'treatment_duration_weeks'. You can override using the `.groups` argument.
# Display table
kbl(descriptive_stats, caption = "Descriptive statistics") %>%
  kable_styling(bootstrap_options = "striped")
Descriptive statistics
treatment_duration_weeks treatment mean sd sem
3 Vehicle 100.000000 0.000000 0.0000000
3 MSH3aso-0.022uM 103.753600 42.175955 24.3502992
3 MSH3aso-0.26uM 19.315167 3.901393 2.2524702
3 MSH3aso-3uM 5.976033 7.204422 4.1594753
3 SCRaso-3uM 190.146533 99.644166 57.5295859
9 Vehicle 100.000000 0.000000 0.0000000
9 MSH3aso-0.022uM 50.742850 39.062771 27.6215500
9 MSH3aso-0.26uM 47.476467 35.675978 20.5975358
9 MSH3aso-3uM 1.235200 NA NA
9 SCRaso-3uM 159.393400 51.249458 29.5888885
12 Vehicle 100.000000 0.000000 0.0000000
12 MSH3aso-0.022uM 34.276467 9.674190 5.5853963
12 MSH3aso-0.26uM 20.043367 23.866658 13.7794212
12 MSH3aso-3uM 0.678400 1.112659 0.6423939
12 SCRaso-3uM 150.230967 36.021422 20.7969779
15 Vehicle 100.000000 0.000000 0.0000000
15 MSH3aso-0.022uM 46.574467 55.495055 32.0400852
15 MSH3aso-0.26uM 23.504833 36.844000 21.2718932
15 MSH3aso-3uM 3.307650 4.641237 3.2818500
15 SCRaso-3uM 132.536300 54.266488 31.3307716

Group mean

This fits a repeated measures linear mixed-effects model ANOVA, lmer(level ~ treatment * treatment_duration_weeks + (1|rep)), to control for changes over time and different biological replicates. It calculates the mean MSH3 levels for each treatment using estimated marginal means and performs posthoc pairwise comparisons. The results are displayed in a table, and the group means are visualised in a plot.

# Repeated measures linear mixed-effects model ANOVA
my_lm <- lmer(level ~ treatment * treatment_duration_weeks + (1|rep), data = plot_data)
anova(my_lm)
# Number of timecourses contributing to each group
sample_count <- plot_data %>%
  distinct(treatment_duration_weeks, treatment) %>%
  group_by(treatment) %>%
  dplyr::summarise(n_timepoints = n()) %>%
  ungroup()

# Posthoc test calculating and comparing group averages
my_emm <- emmeans(my_lm, pairwise ~ treatment | treatment_duration_weeks) # Post-hoc tests of group averages

# Extract group means
group_mean <- as.data.frame(summary(my_emm)$emmeans) %>%
  left_join(sample_count, by = join_by(treatment)) %>%
  dplyr::mutate(SD = SE * sqrt(n_timepoints)) %>%
  relocate(SD, .before = "SE")

# Display group means
kbl(group_mean, caption = "Mean MSH3 level for each treatment group") %>%
  kable_styling(bootstrap_options = "striped")
Mean MSH3 level for each treatment group
treatment treatment_duration_weeks emmean SD SE df lower.CL upper.CL n_timepoints
Vehicle 9.7 100.000000 18.08702 9.043511 14.36387 80.649594 119.35041 4
MSH3aso-0.022uM 9.7 61.291940 21.64865 10.824327 21.37144 38.805318 83.77856 4
MSH3aso-0.26uM 9.7 28.036304 20.63006 10.315032 19.98939 6.518792 49.55382 4
MSH3aso-3uM 9.7 3.497054 23.82746 11.913730 26.05924 -20.989259 27.98337 4
SCRaso-3uM 9.7 158.764153 20.63006 10.315032 19.98939 137.246642 180.28166 4
# Pairwise comparison of treatment groups
group_mean_compare <- as.data.frame(summary(my_emm)$contrasts) %>%
  dplyr::rename(p.adj = p.value) %>% # emmeans already adjusts for multiple comparisons using tukey
  tidyr::separate(contrast, c("group1", "group2"), sep=" - ") %>%
  dplyr::mutate(p.signif = symnum(p.adj,
                                  corr = FALSE,
                                  na = FALSE,
                                  legend = FALSE,
                                  cutpoints = c(0, 0.001, 0.01, 0.05, 0.1, 1),
                                  symbols = c("***", "**", "*", ".", "")),
                group1 = gsub("^\\(|\\)$", "", group1), # Remove brackets
                group2 = gsub("^\\(|\\)$", "", group2)) # Remove brackets

# Hide non-significant results
if (hide_ns) {
  group_mean_compare_plot <- group_mean_compare %>%
    filter(as.numeric(p.adj) < 0.05)
  } else {
    group_mean_compare_plot <- pwc_result
  }

# Display pairwise comparisons
kbl(group_mean_compare, caption = "Pairwise comparison of treatment group means") %>%
  kable_styling(bootstrap_options = "striped")
Pairwise comparison of treatment group means
group1 group2 treatment_duration_weeks estimate SE df t.ratio p.adj p.signif
Vehicle MSH3aso-0.022uM 9.7 38.70806 12.91980 48.38705 2.996025 0.0333704
Vehicle MSH3aso-0.26uM 9.7 71.96370 12.49621 47.78943 5.758844 0.0000057 ***
Vehicle MSH3aso-3uM 9.7 96.50295 13.84530 48.69702 6.970089 0.0000001 ***
Vehicle SCRaso-3uM 9.7 -58.76415 12.49621 47.78943 -4.702560 0.0002079 ***
MSH3aso-0.022uM MSH3aso-0.26uM 9.7 33.25564 13.57150 47.29441 2.450402 0.1198478
MSH3aso-0.022uM MSH3aso-3uM 9.7 57.79489 14.64059 47.28273 3.947579 0.0023225 **
MSH3aso-0.022uM SCRaso-3uM 9.7 -97.47221 13.57150 47.29441 -7.182124 0.0000000 ***
MSH3aso-0.26uM MSH3aso-3uM 9.7 24.53925 14.39221 47.44504 1.705037 0.4407232
MSH3aso-0.26uM SCRaso-3uM 9.7 -130.72785 13.24338 47.09175 -9.871183 0.0000000 ***
MSH3aso-3uM SCRaso-3uM 9.7 -155.26710 14.39221 47.44504 -10.788273 0.0000000 ***
# Plot
my_plot <-
  ggplot(group_mean, aes(x = treatment, y = emmean)) +
  geom_bar(stat = "identity",
           colour = "black",
           # aes(fill = treatment), alpha = 0.8) +
           aes(fill = treatment)) +
  geom_jitter(data = plot_data, 
              aes(x = treatment, y = level, fill = treatment),
              height = 0, width = 0.1,
              size = 3, shape = 21, colour = "black", alpha = 0.8) +
  geom_errorbar(aes(ymin = case_when(errorbar_metric == "CL" ~ lower.CL,
                                     errorbar_metric == "SE" ~ emmean - SE,
                                     errorbar_metric == "SD" ~ emmean - SD,
                                     TRUE ~ NA),
                    ymax = case_when(errorbar_metric == "CL" ~ upper.CL,
                                     errorbar_metric == "SE" ~ emmean + SE,
                                     errorbar_metric == "SD" ~ emmean + SD,
                                     TRUE ~ NA)),
                width = 0.2) +
  stat_pvalue_manual(group_mean_compare_plot,
                     label = "{p.signif}",
                     y.position = ifelse(errorbar_metric == "CL",
                                         1.1 * max(group_mean$upper.CL),
                                         1.1 * (max(group_mean$emmean) + max(group_mean$SE))),
                     step.increase = 0.1,
                     vjust = 0.75,
                     hide.ns = hide_ns,
                     size = 6) +
  scale_x_discrete(labels = treatment_rename,
                   guide = guide_axis(angle = 45)) +
  scale_y_continuous(breaks = pretty_breaks()) +
  scale_colour_manual(values = treatment_colour) +
  scale_fill_manual(values = treatment_colour,) +
  labs(x = "Treatment",
       y = "Relative MSH3 expression (%)",
       colour = "Treatment",
       fill = "Treatment") +
  theme_bw() +
  theme(text = element_text(size = element_text_size),
        legend.position = "none")

# Dispay plot
print(my_plot)

# Export plot
for (output_format in output_formats) {
  ggsave(filename = file.path(out_dir, experiment_id, paste0(experiment_id, "-pwc_group_means.", output_format)),
         plot = my_plot,
         height = 6, width = 6)
}

Group slope

This calculates and compares the slopes of change in MSH3 levels for each treatment. The slopes are estimated using the emtrends function from the emmeans package, which fits a repeated measures linear mixed-effects model ANOVA. Posthoc pairwise comparisons are performed to compare the slopes between treatment groups. The results are displayed in a table, and the slopes are visualised in a plot.

# Posthoc test calculating and comparing group slopes
my_emt <- emtrends(my_lm, ~ treatment, var = "treatment_duration_weeks")

# Extract group slopes
group_slope <- as.data.frame(my_emt)

# Display group slopes
kbl(group_slope, caption = "Slope of change in MSH3 level for each treatment group") %>%
  kable_styling(bootstrap_options = "striped")
Slope of change in MSH3 level for each treatment group
treatment treatment_duration_weeks.trend SE df lower.CL upper.CL
Vehicle 0.0000000 1.827644 47.09175 -3.676557 3.6765570
MSH3aso-0.022uM -5.0192411 2.158652 48.92396 -9.357385 -0.6810973
MSH3aso-0.26uM 0.3708782 2.160993 49.08288 -3.971614 4.7133708
MSH3aso-3uM -0.0562575 2.289665 48.28564 -4.659234 4.5467194
SCRaso-3uM -4.3492758 2.160993 49.08288 -8.691768 -0.0067832
# Pairwise comparisons
my_pwc <- pairs(my_emt)
group_slope_compare <- data.frame(my_pwc) %>%
  dplyr::rename(p.adj = p.value) %>% # emtrends already adjusts for multiple comparisons using tukey
  tidyr::separate(contrast, c("group1", "group2"), sep=" - ") %>%
  dplyr::mutate(p.signif = symnum(p.adj,
                                  corr = FALSE,
                                  na = FALSE,
                                  legend = FALSE,
                                  cutpoints = c(0, 0.001, 0.01, 0.05, 0.1, 1),
                                  symbols = c("***", "**", "*", ".", "")))

# Display pairwise comparisons
kbl(group_slope_compare, caption = "Pairwise comparison of treatment group slopes") %>%
  kable_styling(bootstrap_options = "striped")
Pairwise comparison of treatment group slopes
group1 group2 estimate SE df t.ratio p.adj p.signif
Vehicle (MSH3aso-0.022uM) 5.0192411 2.828438 48.24414 1.7745629 0.3998857
Vehicle (MSH3aso-0.26uM) -0.3708782 2.830225 48.36083 -0.1310419 0.9999310
Vehicle (MSH3aso-3uM) 0.0562575 2.929650 47.84872 0.0192028 1.0000000
Vehicle (SCRaso-3uM) 4.3492758 2.830225 48.36083 1.5367243 0.5442953
(MSH3aso-0.022uM) (MSH3aso-0.26uM) -5.3901193 2.986661 47.09787 -1.8047309 0.3830559
(MSH3aso-0.022uM) (MSH3aso-3uM) -4.9629835 3.102618 47.32952 -1.5996115 0.5050440
(MSH3aso-0.022uM) (SCRaso-3uM) -0.6699653 2.986661 47.09787 -0.2243192 0.9994170
(MSH3aso-0.26uM) (MSH3aso-3uM) 0.4271357 3.100407 47.32645 0.1377676 0.9999157
(MSH3aso-0.26uM) (SCRaso-3uM) 4.7201540 2.984531 47.09175 1.5815398 0.5163374
(MSH3aso-3uM) (SCRaso-3uM) 4.2930183 3.100407 47.32645 1.3846626 0.6402912
# Display group slopes
plot(my_emt)

Group mean at each timepoint

This performs separate ANOVA (aov(MSH3_level ~ treatment)) and pairwise comparisons with posthoc Tukey’s HSD tests for each time point to compare the treatment groups. The results show pairwise comparisons of treatment group means at each specific time point, providing insights into how treatments differ at each stage of the experiment. The results are displayed in separate tables for each time point.

# Vector of time points
my_timepoints <- unique(plot_data$treatment_duration_weeks)
names(my_timepoints) <- my_timepoints

# Perform ANOVA and Tukey's HSD at each time point
group_mean_compare_timepoints <- lapply(my_timepoints, function(my_timepoint,
                                                                plot_data) {
  
  # Filter data
  my_data <- plot_data %>%
    dplyr::filter(treatment_duration_weeks == my_timepoint)
  
  # ANOVA
  my_aov <- aov(level ~ treatment, data = my_data)
  
  # Posthoc tukey's HSD for pairwise comparisons
  my_pwc <- TukeyHSD(my_aov)
  
  # Tidy results
  my_pwc_tidy <- tidy(my_pwc) %>%
    dplyr::mutate(timepoint = my_timepoint,
                  p.signif = symnum(adj.p.value,
                                    corr = FALSE,
                                    na = FALSE,
                                    legend = FALSE,
                                    cutpoints = c(0, 0.001, 0.01, 0.05, 0.1, 1),
                                    symbols = c("***", "**", "*", ".", ""))) %>%
    relocate(timepoint)
  
  # Output
  return(my_pwc_tidy)
  
},
plot_data = plot_data)

# Edit names
names(group_mean_compare_timepoints) <- paste0("week", names(group_mean_compare_timepoints))


# Display tables
for (table_name in names(group_mean_compare_timepoints)) {
  print(kable(group_mean_compare_timepoints[[table_name]],
              caption = paste(table_name, "pairwise comparison of treatment group means")))
}
## 
## 
## Table: week3 pairwise comparison of treatment group means
## 
## | timepoint|term      |contrast                       | null.value|  estimate|   conf.low| conf.high| adj.p.value|p.signif |
## |---------:|:---------|:------------------------------|----------:|---------:|----------:|---------:|-----------:|:--------|
## |         3|treatment |MSH3aso-0.022uM-Vehicle        |          0|   3.75360| -110.53382| 118.04102|   0.9999662|         |
## |         3|treatment |MSH3aso-0.26uM-Vehicle         |          0| -80.68483| -194.97225|  33.60258|   0.2207749|         |
## |         3|treatment |MSH3aso-3uM-Vehicle            |          0| -94.02397| -208.31138|  20.26345|   0.1251589|         |
## |         3|treatment |SCRaso-3uM-Vehicle             |          0|  90.14653|  -24.14088| 204.43395|   0.1482275|         |
## |         3|treatment |MSH3aso-0.26uM-MSH3aso-0.022uM |          0| -84.43843| -206.61682|  37.73995|   0.2365331|         |
## |         3|treatment |MSH3aso-3uM-MSH3aso-0.022uM    |          0| -97.77757| -219.95595|  24.40082|   0.1400028|         |
## |         3|treatment |SCRaso-3uM-MSH3aso-0.022uM     |          0|  86.39293|  -35.78545| 208.57132|   0.2196206|         |
## |         3|treatment |MSH3aso-3uM-MSH3aso-0.26uM     |          0| -13.33913| -135.51752| 108.83925|   0.9961523|         |
## |         3|treatment |SCRaso-3uM-MSH3aso-0.26uM      |          0| 170.83137|   48.65298| 293.00975|   0.0061940|**       |
## |         3|treatment |SCRaso-3uM-MSH3aso-3uM         |          0| 184.17050|   61.99211| 306.34889|   0.0035671|**       |
## 
## 
## Table: week9 pairwise comparison of treatment group means
## 
## | timepoint|term      |contrast                       | null.value|   estimate|     conf.low| conf.high| adj.p.value|p.signif |
## |---------:|:---------|:------------------------------|----------:|----------:|------------:|---------:|-----------:|:--------|
## |         9|treatment |MSH3aso-0.022uM-Vehicle        |          0| -49.257150| -151.4014689|  52.88717|   0.5009701|         |
## |         9|treatment |MSH3aso-0.26uM-Vehicle         |          0| -52.523533| -142.6063552|  37.55929|   0.3391298|         |
## |         9|treatment |MSH3aso-3uM-Vehicle            |          0| -98.764800| -230.6325487|  33.10295|   0.1628924|         |
## |         9|treatment |SCRaso-3uM-Vehicle             |          0|  59.393400|  -30.6894219| 149.47622|   0.2445285|         |
## |         9|treatment |MSH3aso-0.26uM-MSH3aso-0.022uM |          0|  -3.266383| -110.9359493| 104.40318|   0.9999661|         |
## |         9|treatment |MSH3aso-3uM-MSH3aso-0.022uM    |          0| -49.507650| -193.9615311|  94.94623|   0.7603098|         |
## |         9|treatment |SCRaso-3uM-MSH3aso-0.022uM     |          0| 108.650550|    0.9809841| 216.32012|   0.0479058|*        |
## |         9|treatment |MSH3aso-3uM-MSH3aso-0.26uM     |          0| -46.241267| -182.4336919|  89.95116|   0.7659766|         |
## |         9|treatment |SCRaso-3uM-MSH3aso-0.26uM      |          0| 111.916933|   15.6143459| 208.21952|   0.0235927|*        |
## |         9|treatment |SCRaso-3uM-MSH3aso-3uM         |          0| 158.158200|   21.9657748| 294.35063|   0.0236845|*        |
## 
## 
## Table: week12 pairwise comparison of treatment group means
## 
## | timepoint|term      |contrast                       | null.value|  estimate|    conf.low| conf.high| adj.p.value|p.signif |
## |---------:|:---------|:------------------------------|----------:|---------:|-----------:|---------:|-----------:|:--------|
## |        12|treatment |MSH3aso-0.022uM-Vehicle        |          0| -65.72353| -112.375298| -19.07177|   0.0058676|**       |
## |        12|treatment |MSH3aso-0.26uM-Vehicle         |          0| -79.95663| -126.608397| -33.30487|   0.0013056|**       |
## |        12|treatment |MSH3aso-3uM-Vehicle            |          0| -99.32160| -145.973364| -52.66984|   0.0002043|***      |
## |        12|treatment |SCRaso-3uM-Vehicle             |          0|  50.23097|    3.579203|  96.88273|   0.0333263|*        |
## |        12|treatment |MSH3aso-0.26uM-MSH3aso-0.022uM |          0| -14.23310|  -64.105934|  35.63973|   0.8822168|         |
## |        12|treatment |MSH3aso-3uM-MSH3aso-0.022uM    |          0| -33.59807|  -83.470900|  16.27477|   0.2561342|         |
## |        12|treatment |SCRaso-3uM-MSH3aso-0.022uM     |          0| 115.95450|   66.081666| 165.82733|   0.0000917|***      |
## |        12|treatment |MSH3aso-3uM-MSH3aso-0.26uM     |          0| -19.36497|  -69.237800|  30.50787|   0.7213519|         |
## |        12|treatment |SCRaso-3uM-MSH3aso-0.26uM      |          0| 130.18760|   80.314766| 180.06043|   0.0000309|***      |
## |        12|treatment |SCRaso-3uM-MSH3aso-3uM         |          0| 149.55257|   99.679733| 199.42540|   0.0000080|***      |
## 
## 
## Table: week15 pairwise comparison of treatment group means
## 
## | timepoint|term      |contrast                       | null.value|  estimate|    conf.low| conf.high| adj.p.value|p.signif |
## |---------:|:---------|:------------------------------|----------:|---------:|-----------:|---------:|-----------:|:--------|
## |        15|treatment |MSH3aso-0.022uM-Vehicle        |          0| -53.42553| -150.078839|  43.22777|   0.4144185|         |
## |        15|treatment |MSH3aso-0.26uM-Vehicle         |          0| -76.49517| -173.148473|  20.15814|   0.1426786|         |
## |        15|treatment |MSH3aso-3uM-Vehicle            |          0| -96.69235| -206.286898|  12.90220|   0.0909143|.        |
## |        15|treatment |SCRaso-3uM-Vehicle             |          0|  32.53630|  -64.117006| 129.18961|   0.7990783|         |
## |        15|treatment |MSH3aso-0.26uM-MSH3aso-0.022uM |          0| -23.06963| -126.396364|  80.25710|   0.9432480|         |
## |        15|treatment |MSH3aso-3uM-MSH3aso-0.022uM    |          0| -43.26682| -158.789613|  72.25598|   0.7342527|         |
## |        15|treatment |SCRaso-3uM-MSH3aso-0.022uM     |          0|  85.96183|  -17.364897| 189.28856|   0.1169050|         |
## |        15|treatment |MSH3aso-3uM-MSH3aso-0.26uM     |          0| -20.19718| -135.719980|  95.32561|   0.9758222|         |
## |        15|treatment |SCRaso-3uM-MSH3aso-0.26uM      |          0| 109.03147|    5.704736| 212.35820|   0.0377067|*        |
## |        15|treatment |SCRaso-3uM-MSH3aso-3uM         |          0| 129.22865|   13.705853| 244.75145|   0.0272796|*        |

MSH3ko western

# Experiment ID
experiment_id <- names(MSH3aso_data)[[5]]
experiment_id
## [1] "MSH3ko_western_MSH2_MSH6"
# Make results directory
dir.create(file.path(out_dir, experiment_id),
           recursive = TRUE)

# Format data
plot_data <- MSH3aso_data[[experiment_id]] %>%
  pivot_longer(!protein,
               names_to = "label",
               values_to = "level") %>%
  separate_wider_delim(label, "_", names = c("genotype", "rep")) %>%
  dplyr::mutate(level = as.numeric(level) * 100) %>%
  dplyr::filter(!is.na(level))

# Relevel variables
plot_data <- plot_data %>%
  dplyr::mutate(genotype = fct_relevel(genotype, intersect(genotype_order, unique(plot_data$genotype))))

# Comparisons for stat_compare_means
pairwise_comparisons <- lapply(combn(unique(plot_data$genotype), 2, simplify = FALSE), as.character)

# t-test with correction
ttest <- 
  plot_data %>%
  group_by(protein) %>%
  t_test(level ~ genotype) %>%
  adjust_pvalue(method = padj_method) %>%
  add_significance(cutpoints = c(0, 0.001, 0.01, 0.05, 1),
                   symbols = c("***", "**", "*", "ns")) %>%
  add_y_position() %>%
  ungroup()

# Display t-test
kbl(ttest, caption = "T-test") %>%
  kable_styling(bootstrap_options = "striped")
T-test
protein .y. group1 group2 n1 n2 statistic df p p.adj p.adj.signif y.position groups
MSH2 level Unedited MSH3ko 12 9 1.5826844 18.42727 0.131 0.262 ns 176.8654 Unedited, MSH3ko
MSH6 level Unedited MSH3ko 12 9 0.1298626 18.34958 0.898 0.898 ns 224.7324 Unedited, MSH3ko
# ANOVA and corrected posthoc pairwise comparisons
my_aov <- aov(level ~ genotype * protein, data = plot_data)
my_emmeans <- emmeans(my_aov, pairwise ~ genotype | protein) # emmeans 'adjust' argument is giving odd results, with some p values getting more significant
my_pwc <- as.data.frame(my_emmeans$contrasts) %>%
  dplyr::mutate(group1 = str_split_fixed(contrast, " - ", 2)[, 1],
                group2 = str_split_fixed(contrast, " - ", 2)[, 2],
                p.adj = p.adjust(p.value, method = padj_method),
                p.adj.signif = symnum(p.adj, corr = FALSE, na = FALSE, legend = FALSE,
                                      cutpoints = c(0, 0.001, 0.01, 0.05, 1), 
                                      symbols = c("***", "**", "*", "ns"))) %>%
  relocate(c(group1, group2), .after = "contrast") %>%
  relocate(c(p.adj, p.adj.signif), .after = "p.value")
if (hide_ns) {
  my_pwc_plot <- my_pwc %>%
    dplyr::filter(p.adj < 0.05)
} else {
  my_pwc_plot <- my_pwc
}
y_max <- plot_data %>%
  group_by(protein) %>%
  dplyr::summarise(ymax = max(level, na.rm = TRUE))
if (nrow(my_pwc_plot) > 0) {
  my_pwc_plot <- my_pwc_plot %>%
    dplyr::left_join(y_max, by = join_by(protein)) %>%
    group_by(protein) %>%
    dplyr::mutate(y.position = ymax + seq(0.1 * unique(ymax), 0.5 * unique(ymax), length.out = n())) %>%
    ungroup()
} else {
  my_pwc_plot[1,] <- NA
  my_pwc_plot$y.position <- NA
}

# Display ANOVA posthoc
kbl(my_pwc, caption = "ANOVA with corrected posthoc pairwise comparisons") %>%
  kable_styling(bootstrap_options = "striped")
ANOVA with corrected posthoc pairwise comparisons
contrast group1 group2 protein estimate SE df t.ratio p.value p.adj p.adj.signif
Unedited - MSH3ko Unedited MSH3ko MSH2 18.37798 16.42012 38 1.1192356 0.2700647 0.5401294 ns
Unedited - MSH3ko Unedited MSH3ko MSH6 2.50350 16.42012 38 0.1524654 0.8796266 0.8796266 ns
# Plot
my_plot <-
  ggplot(plot_data, aes(x = genotype, y = level)) +
  facet_wrap(vars(protein)) +
  stat_summary(aes(fill = genotype),
               # fun = mean, colour = "black", geom = "bar", alpha = 0.8) +
               fun = mean, colour = "black", geom = "bar") +
  stat_summary(fun.data = mean_se,  
               geom = "errorbar", width = 0.2) +
  geom_jitter(aes(fill = genotype),
              width = 0.2, height = 0, size = point_size, colour = "black", shape = 21) +
  # geom_violin(trim = TRUE) +
  # geom_boxplot(outlier.shape = NA, width = 0.2) +
  # geom_jitter(aes(colour = genotype, fill = genotype),
  #             width = 0.2, height = 0, size = point_size, alpha = 0.75) +
  # stat_summary(fun = mean, geom = "point", shape = 18, size = point_size + 2, show.legend = FALSE) +
  # stat_summary(fun = mean, geom = "text", show.legend = FALSE, vjust = -1, size = point_size + 2,
  #              aes(label = after_stat(round(y, 1)))) +
  # stat_compare_means(comparisons = pairwise_comparisons,
  #                    label = "p.signif",
  #                    hide.ns = hide_ns) +
  # stat_pvalue_manual(ttest, hide.ns = hide_ns, label = "p.adj") +
  stat_pvalue_manual(my_pwc_plot, hide.ns = hide_ns, label = "p.adj.signif", tip.length = 0.01) +
  geom_hline(yintercept = 0, colour = "darkgrey") +
  geom_hline(yintercept = 100, colour = "darkgrey", linetype = "dashed") +
  scale_colour_manual(values = genotype_colour) +
  scale_fill_manual(values = genotype_colour) +
  scale_x_discrete(labels = genotype_rename,
                   guide = guide_axis(angle = 45)) +
  scale_y_continuous(breaks = pretty_breaks()) +
  labs(x = "Genotype",
       y = "Relative protein level (%)",
       colour = "Genotype",
       fill = "Genotype") +
  theme_bw() +
  theme(text = element_text(size = element_text_size),
        legend.position = "none")
my_plot

# Export
for (output_format in output_formats) {
  ggsave(filename = file.path(out_dir, experiment_id, paste0(experiment_id, ".", output_format)),
         plot = my_plot,
         height = 6, width = 6)
}

Small pool PCR

Data from FAN1ko and unedited expansion and contraction tables has been combined into a single facetted plot.

Analysis

# Experiment IDs
experiment_ids <- setNames(names(MSH3aso_data)[6:9], names(MSH3aso_data)[6:9])

# Make results directory
dir.create(file.path(out_dir, "sppcr"),
           recursive = TRUE)

# Format data
plot_data <- lapply(experiment_ids, function(experiment_id,
                                             MSH3aso_data,
                                             treatment_rename) {
   
  # Extract annotation from experiment ID
  experiment_id_components <- as.character(str_split(experiment_id, "_", simplify = TRUE))
  my_genotype = experiment_id_components[[2]]
  my_metric = experiment_id_components[[3]]
  
  # Format
  my_data <- MSH3aso_data[[experiment_id]] %>%
    pivot_longer(!Metric,
               names_to = "label",
               values_to = "level") %>%
    separate_wider_delim(label, "_", names = c("treatment", "rep")) %>%
    dplyr::mutate(level = as.numeric(level)) %>%
    dplyr::filter(!is.na(level)) %>%
    dplyr::mutate(genotype = my_genotype,
                  metric = my_metric,
                  genotype = ifelse(genotype == "unedited", "Unedited", genotype)) %>%
    relocate(c(metric, genotype), .after = "Metric")
  
  # Output
  return(my_data)
  
},
MSH3aso_data = MSH3aso_data,
treatment_rename = treatment_rename)

# rbind
plot_data <- rbindlist(plot_data)

# Exclude SCR treatment from unedited genotype
plot_data <- plot_data %>%
  dplyr::filter(!(genotype == "Unedited" & treatment == "SCRaso-3uM"))

# Relevel
plot_data <- plot_data %>%
  dplyr::mutate(treatment = fct_relevel(treatment, intersect(treatment_order, unique(plot_data$treatment))),
                # genotype = fct_relevel(genotype, intersect(genotype_order, unique(plot_data$genotype))),
                genotype = fct_relevel(genotype, c("FAN1ko", "Unedited")),
                Metric = fct_relevel(Metric, c("Expansions per lane", "Contractions per lane")))

# Add to treatment_rename
my_treatment_rename <- sapply(names(treatment_rename), function(x) {
  if (x == "Baseline") {
    paste(treatment_rename[x], "(36d)")
  } else {
    paste(treatment_rename[x], "(12wks)")
  }
}, USE.NAMES = TRUE)

# t-test
ttest <- 
  plot_data %>%
  group_by(Metric, genotype) %>%
  t_test(level ~ treatment) %>%
  adjust_pvalue(method = padj_method) %>%
  add_significance(cutpoints = c(0, 0.001, 0.01, 0.05, 1),
                   symbols = c("***", "**", "*", "ns")) %>%
  add_y_position() %>%
  ungroup()

# Display t-test
kbl(ttest, caption = "T-test") %>%
  kable_styling(bootstrap_options = "striped")
T-test
Metric genotype .y. group1 group2 n1 n2 statistic df p p.adj p.adj.signif y.position groups
Expansions per lane FAN1ko level Baseline Vehicle 3 3 -1.6223360 3.145469 0.199 0.8516250 ns 0.9830 Baseline, Vehicle
Expansions per lane FAN1ko level Baseline MSH3aso-3uM 3 3 0.3630952 3.997965 0.735 0.8516250 ns 1.1126 Baseline , MSH3aso-3uM
Expansions per lane FAN1ko level Baseline SCRaso-3uM 3 3 -1.3874776 3.032716 0.258 0.8516250 ns 1.2422 Baseline , SCRaso-3uM
Expansions per lane FAN1ko level Vehicle MSH3aso-3uM 3 3 2.1007572 3.188233 0.121 0.8516250 ns 1.3718 Vehicle , MSH3aso-3uM
Expansions per lane FAN1ko level Vehicle SCRaso-3uM 3 3 0.3786803 3.984831 0.724 0.8516250 ns 1.5014 Vehicle , SCRaso-3uM
Expansions per lane FAN1ko level MSH3aso-3uM SCRaso-3uM 3 3 -1.8691854 3.073112 0.156 0.8516250 ns 1.6310 MSH3aso-3uM, SCRaso-3uM
Expansions per lane Unedited level Baseline Vehicle 3 3 -0.5341094 3.172790 0.628 0.8516250 ns 1.1080 Baseline, Vehicle
Expansions per lane Unedited level Baseline MSH3aso-3uM 3 3 -0.0960038 2.991441 0.930 0.9847059 ns 1.2700 Baseline , MSH3aso-3uM
Expansions per lane Unedited level Vehicle MSH3aso-3uM 3 3 0.3322949 3.960746 0.757 0.8516250 ns 1.4320 Vehicle , MSH3aso-3uM
Contractions per lane FAN1ko level Baseline Vehicle 3 3 -0.6167433 3.079932 0.580 0.8516250 ns 0.6790 Baseline, Vehicle
Contractions per lane FAN1ko level Baseline MSH3aso-3uM 3 3 0.4007954 3.967370 0.709 0.8516250 ns 0.8086 Baseline , MSH3aso-3uM
Contractions per lane FAN1ko level Baseline SCRaso-3uM 3 3 0.0194998 3.999821 0.985 0.9850000 ns 0.9382 Baseline , SCRaso-3uM
Contractions per lane FAN1ko level Vehicle MSH3aso-3uM 3 3 0.8921462 2.922612 0.440 0.8516250 ns 1.0678 Vehicle , MSH3aso-3uM
Contractions per lane FAN1ko level Vehicle SCRaso-3uM 3 3 0.6289575 3.092149 0.573 0.8516250 ns 1.1974 Vehicle , SCRaso-3uM
Contractions per lane FAN1ko level MSH3aso-3uM SCRaso-3uM 3 3 -0.3789710 3.962464 0.724 0.8516250 ns 1.3270 MSH3aso-3uM, SCRaso-3uM
Contractions per lane Unedited level Baseline Vehicle 3 3 0.9637981 2.789449 0.411 0.8516250 ns 0.9410 Baseline, Vehicle
Contractions per lane Unedited level Baseline MSH3aso-3uM 3 3 -0.5800253 2.278657 0.614 0.8516250 ns 1.1030 Baseline , MSH3aso-3uM
Contractions per lane Unedited level Vehicle MSH3aso-3uM 3 3 -0.8737264 2.057593 0.472 0.8516250 ns 1.2650 Vehicle , MSH3aso-3uM
# ANOVA and corrected posthoc pairwise comparisons
my_aov <- aov(level ~ genotype * Metric, data = plot_data)
my_emmeans <- emmeans(my_aov, pairwise ~ genotype | Metric) # emmeans 'adjust' argument is giving odd results, with some p values getting more significant
my_pwc <- as.data.frame(my_emmeans$contrasts) %>%
  dplyr::mutate(group1 = str_split_fixed(contrast, " - ", 2)[, 1],
                group2 = str_split_fixed(contrast, " - ", 2)[, 2],
                p.adj = p.adjust(p.value, method = padj_method),
                p.adj.signif = symnum(p.adj, corr = FALSE, na = FALSE, legend = FALSE,
                                      cutpoints = c(0, 0.001, 0.01, 0.05, 1), 
                                      symbols = c("***", "**", "*", "ns"))) %>%
  relocate(c(group1, group2), .after = "contrast") %>%
  relocate(c(p.adj, p.adj.signif), .after = "p.value")
if (hide_ns) {
  my_pwc_plot <- my_pwc %>%
    dplyr::filter(p.adj < 0.05)
} else {
  my_pwc_plot <- my_pwc
}
y_max <- plot_data %>%
  group_by(Metric) %>%
  dplyr::summarise(ymax = max(level, na.rm = TRUE))
if (nrow(my_pwc_plot) > 0) {
  my_pwc_plot <- my_pwc_plot %>%
    dplyr::left_join(y_max, by = join_by(Metric)) %>%
    group_by(Metric) %>%
    dplyr::mutate(y.position = ymax + seq(0.1 * unique(ymax), 0.5 * unique(ymax), length.out = n())) %>%
    ungroup()
} else {
  my_pwc_plot[1,] <- NA
  my_pwc_plot$y.position <- NA
}

# Display ANOVA posthoc
kbl(my_pwc, caption = "ANOVA with corrected posthoc pairwise comparisons") %>%
  kable_styling(bootstrap_options = "striped")
ANOVA with corrected posthoc pairwise comparisons
contrast group1 group2 Metric estimate SE df t.ratio p.value p.adj p.adj.signif
FAN1ko - Unedited FAN1ko Unedited Expansions per lane 0.0922244 0.109188 38 0.8446387 0.4035991 0.8071983 ns
FAN1ko - Unedited FAN1ko Unedited Contractions per lane 0.0086691 0.109188 38 0.0793959 0.9371344 0.9371344 ns

Plot together

# Plot
my_plot <-
  ggplot(plot_data, aes(x = treatment, y = level)) +
  facet_grid(cols = vars(Metric),
             rows = vars(genotype),
             labeller = labeller(genotype = genotype_rename)) +
  geom_hline(yintercept = 0, colour = "darkgrey") +
  stat_summary(aes(fill = treatment),
               # fun = mean, colour = "black", geom = "bar", alpha = 0.8) +
               fun = mean, colour = "black", geom = "bar") +
  stat_summary(fun.data = mean_se,  
               geom = "errorbar", width = 0.2) +
  geom_jitter(aes(fill = treatment),
              width = 0.2, height = 0, size = point_size, colour = "black", shape = 21) +
  # geom_violin(trim = TRUE) +
  # geom_boxplot(outlier.shape = NA, width = 0.2) +
  # geom_jitter(aes(colour = treatment, fill = treatment),
  #             width = 0.1, height = 0, size = point_size, alpha = 0.8) +
  # stat_summary(fun = mean, geom = "point", shape = 18, size = point_size + 2, show.legend = FALSE) +
  # stat_summary(fun = mean, geom = "text", show.legend = FALSE, vjust = -1, size = point_size + 2,
  #              aes(label = after_stat(round(y, 1)))) +
  stat_pvalue_manual(ttest, hide.ns = hide_ns, label = "p.adj", tip.length = 0.01) +
  scale_colour_manual(values = treatment_colour) +
  scale_fill_manual(values = treatment_colour) +
  scale_x_discrete(labels = my_treatment_rename,
                   guide = guide_axis(angle = 45)) +
  scale_y_continuous(breaks = pretty_breaks()) +
  labs(x = "Treatment",
       y = "Instability events per lane",
       colour = "Treatment",
       fill = "Treatment") +
  theme_bw() +
  theme(text = element_text(size = element_text_size),
        legend.position = "none")
my_plot

# Export
for (output_format in output_formats) {
  ggsave(filename = file.path(out_dir, "sppcr", paste0(experiment_id, ".", output_format)),
         plot = my_plot,
         height = 12, width = 12)
}

Plot genotypes separately

# Genotype vector
my_genotypes <- as.character(unique(plot_data$genotype))
names(my_genotypes) <- my_genotypes

# Plot
my_plots <- lapply(my_genotypes, function(my_genotype,
                                          plot_data,
                                          ttest) {
  # Filter data
  my_data <- plot_data %>%
    dplyr::filter(genotype == my_genotype)
  
  # Plot
  my_plot <-
    ggplot(my_data, aes(x = treatment, y = level)) +
    facet_wrap(vars(Metric)) +
    geom_hline(yintercept = 0, colour = "darkgrey") +
    stat_summary(aes(fill = treatment),
                 # fun = mean, colour = "black", geom = "bar", alpha = 0.8) +
                 fun = mean, colour = "black", geom = "bar") +
    stat_summary(fun.data = mean_se,  
                 geom = "errorbar", width = 0.2) +
    geom_jitter(aes(fill = treatment),
                width = 0.2, height = 0, size = point_size, colour = "black", shape = 21) +
    stat_pvalue_manual(ttest, hide.ns = hide_ns, label = "p.adj", tip.length = 0.01) +
    scale_colour_manual(values = treatment_colour) +
    scale_fill_manual(values = treatment_colour) +
    scale_x_discrete(labels = my_treatment_rename,
                     guide = guide_axis(angle = 45)) +
    scale_y_continuous(breaks = pretty_breaks()) +
    labs(x = "Treatment",
         y = "Instability events per lane",
         colour = "Treatment",
         fill = "Treatment") +
    theme_bw() +
    theme(text = element_text(size = element_text_size),
          legend.position = "none")
  
  # Output
  return(my_plot)
  
},
plot_data = plot_data,
ttest = ttest)

# Display plots
for (my_plot in my_plots) {
  print(my_plot)
}

# Export
for (plot_name in names(my_plots)) {
  my_plot <- my_plots[[plot_name]]
  for (output_format in output_formats) {
    ggsave(filename = file.path(out_dir, "sppcr", paste0(experiment_id, paste0("_", plot_name, "."), output_format)),
           plot = my_plot,
           height = 6, width = 10)
  }
}

Plot genotypes and metrics separately

# Genotype vector
my_genotypes <- as.character(unique(plot_data$genotype))
names(my_genotypes) <- my_genotypes

# Metric vector
my_metrics <- as.character(unique(plot_data$Metric))
names(my_metrics) <- my_metrics

# Plot
my_plots <- lapply(my_genotypes, function(my_genotype,
                                          plot_data,
                                          ttest,
                                          my_metrics) {
  
  # Plot for each metric
  my_plot_metrics <- lapply(my_metrics, function(my_metric,
                                                 my_genotype) {
    
    # # Manual
    # my_metric = my_metrics[[1]]
    
    # Filter data
    my_data <- plot_data %>%
      dplyr::filter(genotype == my_genotype,
                    Metric == my_metric) 
    
    # Plot
    my_plot <-
      ggplot(my_data, aes(x = treatment, y = level)) +
      geom_hline(yintercept = 0, colour = "darkgrey") +
      stat_summary(aes(fill = treatment),
                   # fun = mean, colour = "black", geom = "bar", alpha = 0.8) +
                   fun = mean, colour = "black", geom = "bar") +
      stat_summary(fun.data = mean_se,  
                   geom = "errorbar", width = 0.2) +
      geom_jitter(aes(fill = treatment),
                  width = 0.2, height = 0, size = point_size, colour = "black", shape = 21) +
      stat_pvalue_manual(ttest, hide.ns = hide_ns, label = "p.adj", tip.length = 0.01) +
      scale_colour_manual(values = treatment_colour) +
      scale_fill_manual(values = treatment_colour) +
      scale_x_discrete(labels = my_treatment_rename,
                       guide = guide_axis(angle = 45)) +
      scale_y_continuous(breaks = pretty_breaks()) +
      labs(x = "Treatment",
           y = "Instability events per lane",
           colour = "Treatment",
           fill = "Treatment") +
      theme_bw() +
      theme(text = element_text(size = element_text_size),
            legend.position = "none")
    
    # Output
    return(my_plot)
    
  },
  my_genotype = my_genotype)
  
  # Output
  return(my_plot_metrics)
  
},
plot_data = plot_data,
ttest = ttest,
my_metrics = my_metrics)


# Flatten
my_plots <- unlist(my_plots, recursive = FALSE)

# Display plots
for (my_plot in my_plots) {
  print(my_plot)
}

# Export
for (plot_name in names(my_plots)) {
  my_plot <- my_plots[[plot_name]]
  for (output_format in output_formats) {
    ggsave(filename = file.path(out_dir, "sppcr", paste0(experiment_id, paste0("_", plot_name, "."), output_format)),
           plot = my_plot,
           height = 6, width = 6)
  }
}

In vivo titration

Import data

# Experiment IDs
experiment_ids <- setNames(names(MSH3aso_data)[10:11], names(MSH3aso_data)[10:11])

# Make results directory
dir.create(file.path(out_dir, "in_vivo_titration"),
           recursive = TRUE)

# Format data
plot_data <- lapply(experiment_ids, function(experiment_id,
                                             MSH3aso_data) {
  
  # Format
  my_data <- MSH3aso_data[[experiment_id]] %>%
    pivot_longer(!tissue,
               names_to = "label",
               values_to = "level") %>%
    separate_wider_delim(label, "_", names = c("treatment", "rep")) %>%
    dplyr::mutate(level = as.numeric(level),
                  experiment_id = experiment_id) %>%
    dplyr::filter(!is.na(level)) %>%
    relocate(experiment_id)
    
  # Output
  return(my_data)
  
},
MSH3aso_data = MSH3aso_data)

# rbind
plot_data <- rbindlist(plot_data)

# Relevel
plot_data <- plot_data %>%
  dplyr::mutate(tissue = fct_relevel(tissue, intersect(tissue_order, unique(plot_data$tissue))),
                treatment = fct_relevel(treatment, intersect(treatment, unique(plot_data$treatment))))

Experiments combined

Data from two experiments have been combined for stats and plotting.

# t-test
ttest <- 
  plot_data %>%
  group_by(tissue) %>%
  t_test(level ~ treatment) %>%
  adjust_pvalue(method = padj_method) %>%
  add_significance(cutpoints = c(0, 0.001, 0.01, 0.05, 1),
                   symbols = c("***", "**", "*", "ns")) %>%
  add_y_position() %>%
  ungroup()

# Display t-test
kbl(ttest, caption = "T-test") %>%
  kable_styling(bootstrap_options = "striped")
T-test
tissue .y. group1 group2 n1 n2 statistic df p p.adj p.adj.signif y.position groups
Cortex level Vehicle 3ug 8 8 1.5415216 9.132123 1.57e-01 0.1860741 ns 119.3469 Vehicle, 3ug
Cortex level Vehicle 10ug 8 8 2.5865314 8.569527 3.00e-02 0.0417391
130.4767 Vehicle, 10ug
Cortex level Vehicle 30ug 8 8 13.0916812 9.896674 1.00e-07 0.0000023 *** 141.6066 Vehicle, 30ug
Cortex level Vehicle SCRaso 8 4 -0.8948261 5.464449 4.09e-01 0.4362667 ns 152.7365 Vehicle, SCRaso
Cortex level 3ug 10ug 8 8 1.0006693 13.661257 3.34e-01 0.3685517 ns 163.8663 3ug , 10ug
Cortex level 3ug 30ug 8 8 8.0797152 13.639093 1.50e-06 0.0000088 *** 174.9962 3ug , 30ug
Cortex level 3ug SCRaso 8 4 -1.9654142 9.987477 7.80e-02 0.0998400 ns 186.1261 3ug , SCRaso
Cortex level 10ug 30ug 8 8 6.2577056 12.757057 3.20e-05 0.0000853 *** 197.2559 10ug, 30ug
Cortex level 10ug SCRaso 8 4 -2.9044349 9.731249 1.60e-02 0.0232727
208.3858 10ug , SCRaso
Cortex level 30ug SCRaso 8 4 -12.2273936 9.858507 3.00e-07 0.0000030 *** 219.5157 30ug , SCRaso
Striatum level Vehicle 3ug 4 4 0.5318878 4.900537 6.18e-01 0.6180000 ns 120.4869 Vehicle, 3ug
Striatum level Vehicle 10ug 4 4 2.2498355 3.606236 9.50e-02 0.1169231 ns 132.5071 Vehicle, 10ug
Striatum level Vehicle 30ug 4 5 7.1738608 4.011652 2.00e-03 0.0035556 ** 144.5274 Vehicle, 30ug
Striatum level 3ug 10ug 4 4 2.5708958 4.585510 5.40e-02 0.0720000 ns 156.5476 3ug , 10ug
Striatum level 3ug 30ug 4 5 9.8404489 5.570527 9.93e-05 0.0002444 *** 168.5679 3ug , 30ug
Striatum level 10ug 30ug 4 5 10.3652603 6.925037 1.82e-05 0.0000529 *** 180.5882 10ug, 30ug
Brainstem level Vehicle 3ug 4 4 6.2762166 5.062854 1.00e-03 0.0018824 ** 113.6829 Vehicle, 3ug
Brainstem level Vehicle 10ug 4 4 17.6634697 5.767722 3.00e-06 0.0000122 *** 125.7031 Vehicle, 10ug
Brainstem level Vehicle 30ug 4 5 32.4960641 6.975689 0.00e+00 0.0000002 *** 137.7234 Vehicle, 30ug
Brainstem level 3ug 10ug 4 4 8.0782464 5.643734 2.62e-04 0.0005989 *** 149.7436 3ug , 10ug
Brainstem level 3ug 30ug 4 5 18.8693737 5.569056 2.90e-06 0.0000122 *** 161.7639 3ug , 30ug
Brainstem level 10ug 30ug 4 5 12.2559102 6.544504 9.40e-06 0.0000299 *** 173.7842 10ug, 30ug
Spinal cord level Vehicle 3ug 8 8 3.1598645 9.857447 1.00e-02 0.0152381
117.3129 Vehicle, 3ug
Spinal cord level Vehicle 10ug 8 8 5.6296674 7.974333 4.99e-04 0.0009980 *** 128.4427 Vehicle, 10ug
Spinal cord level Vehicle 30ug 8 8 12.3280829 8.050690 1.70e-06 0.0000088 *** 139.5726 Vehicle, 30ug
Spinal cord level Vehicle SCRaso 8 4 1.4272505 3.660822 2.33e-01 0.2662857 ns 150.7025 Vehicle, SCRaso
Spinal cord level 3ug 10ug 8 8 3.3243802 11.143146 7.00e-03 0.0112000
161.8323 3ug , 10ug
Spinal cord level 3ug 30ug 8 8 9.2092605 11.401101 1.30e-06 0.0000088 *** 172.9622 3ug , 30ug
Spinal cord level 3ug SCRaso 8 4 -0.7888956 6.115341 4.60e-01 0.4748387 ns 184.0921 3ug , SCRaso
Spinal cord level 10ug 30ug 8 8 4.6709975 13.979700 3.62e-04 0.0007723 *** 195.2219 10ug, 30ug
Spinal cord level 10ug SCRaso 8 4 -3.5884061 9.592936 5.00e-03 0.0084211 ** 206.3518 10ug , SCRaso
Spinal cord level 30ug SCRaso 8 4 -8.8324037 9.429196 7.30e-06 0.0000260 *** 217.4817 30ug , SCRaso
# ANOVA and corrected posthoc pairwise comparisons
my_aov <- aov(level ~ treatment * tissue, data = plot_data)
my_emmeans <- emmeans(my_aov, pairwise ~ treatment | tissue) # emmeans 'adjust' argument is giving odd results, with some p values getting more significant
my_pwc <- as.data.frame(my_emmeans$contrasts) %>%
  dplyr::mutate(group1 = str_split_fixed(contrast, " - ", 2)[, 1],
                group2 = str_split_fixed(contrast, " - ", 2)[, 2],
                p.adj = p.adjust(p.value, method = padj_method),
                p.adj.signif = symnum(p.adj, corr = FALSE, na = FALSE, legend = FALSE,
                                      cutpoints = c(0, 0.001, 0.01, 0.05, 1), 
                                      symbols = c("***", "**", "*", "ns"))) %>%
  relocate(c(group1, group2), .after = "contrast") %>%
  relocate(c(p.adj, p.adj.signif), .after = "p.value")
if (hide_ns) {
  my_pwc_plot <- my_pwc %>%
    dplyr::filter(p.adj < 0.05)
} else {
  my_pwc_plot <- my_pwc
}
y_max <- plot_data %>%
  group_by(tissue) %>%
  dplyr::summarise(ymax = max(level, na.rm = TRUE))
if (nrow(my_pwc_plot) > 0) {
  my_pwc_plot <- my_pwc_plot %>%
    dplyr::left_join(y_max, by = join_by(tissue)) %>%
    group_by(tissue) %>%
    dplyr::mutate(y.position = ymax + seq(0.1 * unique(ymax), 0.5 * unique(ymax), length.out = n())) %>%
    ungroup()
} else {
  my_pwc_plot[1,] <- NA
  my_pwc_plot$y.position <- NA
}

# Display ANOVA posthoc
kbl(my_pwc, caption = "ANOVA with corrected posthoc pairwise comparisons") %>%
  kable_styling(bootstrap_options = "striped")
ANOVA with corrected posthoc pairwise comparisons
contrast group1 group2 tissue estimate SE df t.ratio p.value p.adj p.adj.signif
Vehicle - 3ug Vehicle 3ug Cortex 6.688529 5.038409 88 1.3275080 0.6749243 0.7999102 ns
Vehicle - 10ug Vehicle 10ug Cortex 12.910300 5.038409 88 2.5623763 0.0864607 0.1317497 ns
Vehicle - 30ug Vehicle 30ug Cortex 49.453178 5.038409 88 9.8152365 0.0000000 0.0000000 ***
Vehicle - SCRaso Vehicle SCRaso Cortex -2.683435 6.170766 88 -0.4348625 0.9924383 0.9924383 ns
3ug - 10ug 3ug 10ug Cortex 6.221771 5.038409 88 1.2348682 0.7310900 0.8355315 ns
3ug - 30ug 3ug 30ug Cortex 42.764649 5.038409 88 8.4877285 0.0000000 0.0000000 ***
3ug - SCRaso 3ug SCRaso Cortex -9.371964 6.170766 88 -1.5187683 0.5532082 0.6808716 ns
10ug - 30ug 10ug 30ug Cortex 36.542878 5.038409 88 7.2528603 0.0000000 0.0000000 ***
10ug - SCRaso 10ug SCRaso Cortex -15.593735 6.170766 88 -2.5270340 0.0938909 0.1365686 ns
30ug - SCRaso 30ug SCRaso Cortex -52.136613 6.170766 88 -8.4489696 0.0000000 0.0000000 ***
Vehicle - 3ug Vehicle 3ug Striatum 3.659846 7.125387 88 0.5136347 0.9556332 0.9864601 ns
Vehicle - 10ug Vehicle 10ug Striatum 13.950394 7.125387 88 1.9578438 0.2119660 0.2826214 ns
Vehicle - 30ug Vehicle 30ug Striatum 45.805947 6.759735 88 6.7762930 0.0000000 0.0000000 ***
Vehicle - SCRaso Vehicle SCRaso Striatum NA NA NA NA NA NA
3ug - 10ug 3ug 10ug Striatum 10.290548 7.125387 88 1.4442091 0.4755368 0.6086871 ns
3ug - 30ug 3ug 30ug Striatum 42.146101 6.759735 88 6.2348745 0.0000001 0.0000002 ***
3ug - SCRaso 3ug SCRaso Striatum NA NA NA NA NA NA
10ug - 30ug 10ug 30ug Striatum 31.855553 6.759735 88 4.7125445 0.0000532 0.0001064 ***
10ug - SCRaso 10ug SCRaso Striatum NA NA NA NA NA NA
30ug - SCRaso 30ug SCRaso Striatum NA NA NA NA NA NA
Vehicle - 3ug Vehicle 3ug Brainstem 18.896578 7.125387 88 2.6520074 0.0459968 0.0735949 ns
Vehicle - 10ug Vehicle 10ug Brainstem 44.900223 7.125387 88 6.3014438 0.0000001 0.0000002 ***
Vehicle - 30ug Vehicle 30ug Brainstem 77.182277 6.759735 88 11.4179438 0.0000000 0.0000000 ***
Vehicle - SCRaso Vehicle SCRaso Brainstem NA NA NA NA NA NA
3ug - 10ug 3ug 10ug Brainstem 26.003645 7.125387 88 3.6494364 0.0024713 0.0041622 **
3ug - 30ug 3ug 30ug Brainstem 58.285700 6.759735 88 8.6224826 0.0000000 0.0000000 ***
3ug - SCRaso 3ug SCRaso Brainstem NA NA NA NA NA NA
10ug - 30ug 10ug 30ug Brainstem 32.282054 6.759735 88 4.7756388 0.0000416 0.0000887 ***
10ug - SCRaso 10ug SCRaso Brainstem NA NA NA NA NA NA
30ug - SCRaso 30ug SCRaso Brainstem NA NA NA NA NA NA
Vehicle - 3ug Vehicle 3ug Spinal cord 11.934323 5.038409 88 2.3686688 0.1337486 0.1860850 ns
Vehicle - 10ug Vehicle 10ug Spinal cord 34.876693 5.038409 88 6.9221637 0.0000000 0.0000000 ***
Vehicle - 30ug Vehicle 30ug Spinal cord 73.707582 5.038409 88 14.6291378 0.0000000 0.0000000 ***
Vehicle - SCRaso Vehicle SCRaso Spinal cord 7.259351 6.170766 88 1.1764101 0.7647539 0.8438664 ns
3ug - 10ug 3ug 10ug Spinal cord 22.942370 5.038409 88 4.5534949 0.0001613 0.0003037 ***
3ug - 30ug 3ug 30ug Spinal cord 61.773259 5.038409 88 12.2604689 0.0000000 0.0000000 ***
3ug - SCRaso 3ug SCRaso Spinal cord -4.674972 6.170766 88 -0.7576000 0.9419017 0.9864601 ns
10ug - 30ug 10ug 30ug Spinal cord 38.830889 5.038409 88 7.7069741 0.0000000 0.0000000 ***
10ug - SCRaso 10ug SCRaso Spinal cord -27.617342 6.170766 88 -4.4755130 0.0002166 0.0003851 ***
30ug - SCRaso 30ug SCRaso Spinal cord -66.448231 6.170766 88 -10.7682309 0.0000000 0.0000000 ***
# Shape rename
my_shape_rename <- setNames(c("1", "2"), unique(plot_data$experiment_id))

# Plot
my_plot <-
  ggplot(plot_data, aes(x = treatment, y = level)) +
  facet_wrap(vars(tissue), nrow = 1, scales = "free_x") +
  stat_summary(aes(fill = treatment),
               # fun = mean, colour = "black", geom = "bar", alpha = 0.8) +
               fun = mean, colour = "black", geom = "bar") +
  stat_summary(fun.data = mean_se,  
               geom = "errorbar", width = 0.2) +
  geom_jitter(aes(fill = treatment, shape = experiment_id),
              width = 0.2, height = 0, size = point_size, colour = "black") +
  # geom_violin(trim = TRUE) +
  # geom_boxplot(outlier.shape = NA, width = 0.2) +
  # geom_jitter(aes(colour = treatment, fill = treatment, shape = experiment_id),
  #             width = 0.1, height = 0, size = point_size, alpha = 0.8) +
  # stat_summary(fun = mean, geom = "point", shape = 18, size = point_size + 2, show.legend = FALSE) +
  # stat_summary(fun = mean, geom = "text", show.legend = FALSE, vjust = -1, size = point_size + 2,
  #              aes(label = after_stat(round(y, 1)))) +
  geom_hline(yintercept = 0, colour = "darkgrey") +
  geom_hline(yintercept = 100, colour = "darkgrey", linetype = "dashed") +
  # stat_pvalue_manual(ttest, hide.ns = hide_ns, label = "p.adj.signif")
  stat_pvalue_manual(my_pwc_plot, label = "p.adj.signif", hide.ns = TRUE, tip.length = 0.01) +
  scale_x_discrete(labels = treatment_ionis_rename,
                   guide = guide_axis(angle = 45)) +
  scale_y_continuous(breaks = pretty_breaks()) +
  scale_colour_manual(values = treatment_ionis_colour,
                      labels = treatment_ionis_rename) +
  scale_fill_manual(values = treatment_ionis_colour,
                    labels = treatment_ionis_rename) +
  # scale_shape_discrete(labels = my_shape_rename) +
  scale_shape_manual(values = c(21, 24),
                     labels = my_shape_rename) +
  labs(x = "Treatment",
       y = "Relative expression (%)",
       colour = "Treatment",
       fill = "Treatment",
       shape = "Experiment") +
  theme_minimal() +
  theme(text = element_text(size = element_text_size))
my_plot

# Export
for (output_format in output_formats) {
  ggsave(filename = file.path(out_dir, "in_vivo_titration", paste0("in_vivo_titration-combined", ".", output_format)),
         plot = my_plot,
         height = 7, width = 14)
}

Experiment separated

Data from the two experiments kept separate for stats and plotting.

# Vector of experiments
my_experiment_names <- unique(plot_data$experiment_id)
names(my_experiment_names) <- my_experiment_names

# Plot
my_plots <- lapply(my_experiment_names, function(my_experiment_name,
                                                 plot_data,
                                                 padj_method,
                                                 hide_ns,
                                                 treatment_ionis_rename,
                                                 treatment_ionis_colour,
                                                 element_text_size) {
  
  # Filter data
  my_data <- plot_data %>%
    dplyr::filter(experiment_id == my_experiment_name)
  
  # t-test
  ttest <- 
    my_data %>%
    group_by(tissue) %>%
    t_test(level ~ treatment) %>%
    adjust_pvalue(method = padj_method) %>%
    add_significance(cutpoints = c(0, 0.001, 0.01, 0.05, 1),
                     symbols = c("***", "**", "*", "ns")) %>%
    add_y_position() %>%
    ungroup()
  
  
  # ANOVA and corrected posthoc pairwise comparisons
  my_aov <- aov(level ~ treatment * tissue, data = my_data)
  my_emmeans <- emmeans(my_aov, pairwise ~ treatment | tissue) # emmeans 'adjust' argument is giving odd results, with some p values getting more significant
  my_pwc <- as.data.frame(my_emmeans$contrasts) %>%
    dplyr::mutate(group1 = str_split_fixed(contrast, " - ", 2)[, 1],
                  group2 = str_split_fixed(contrast, " - ", 2)[, 2],
                  p.adj = p.adjust(p.value, method = padj_method),
                  p.adj.signif = symnum(p.adj, corr = FALSE, na = FALSE, legend = FALSE,
                                        cutpoints = c(0, 0.001, 0.01, 0.05, 1), 
                                        symbols = c("***", "**", "*", "ns"))) %>%
    relocate(c(group1, group2), .after = "contrast") %>%
    relocate(c(p.adj, p.adj.signif), .after = "p.value")
  if (hide_ns) {
    my_pwc_plot <- my_pwc %>%
      dplyr::filter(p.adj < 0.05)
  } else {
    my_pwc_plot <- my_pwc
  }
  y_max <- my_data %>%
    group_by(tissue) %>%
    dplyr::summarise(ymax = max(level, na.rm = TRUE))
  if (nrow(my_pwc_plot) > 0) {
    my_pwc_plot <- my_pwc_plot %>%
      dplyr::left_join(y_max, by = join_by(tissue)) %>%
      group_by(tissue) %>%
      dplyr::mutate(y.position = ymax + seq(0.1 * unique(ymax), 0.5 * unique(ymax), length.out = n())) %>%
      ungroup()
  } else {
    my_pwc_plot[1,] <- NA
    my_pwc_plot$y.position <- NA
  }
  
  # Plot
  my_plot <-
    ggplot(my_data, aes(x = treatment, y = level)) +
    facet_wrap(vars(tissue), nrow = 1, scales = "free_x") +
    stat_summary(aes(fill = treatment),
                 # fun = mean, colour = "black", geom = "bar", alpha = 0.8) +
                 fun = mean, colour = "black", geom = "bar") +
    stat_summary(fun.data = mean_se,  
                 geom = "errorbar", width = 0.2) +
    geom_jitter(aes(fill = treatment),
                width = 0.2, height = 0, size = point_size, colour = "black", shape = 21) +
    geom_hline(yintercept = 0, colour = "darkgrey") +
    geom_hline(yintercept = 100, colour = "darkgrey", linetype = "dashed") +
    # stat_pvalue_manual(ttest, hide.ns = hide_ns, label = "p.adj.signif")
    stat_pvalue_manual(my_pwc_plot, label = "p.adj.signif", hide.ns = hide_ns, tip.length = 0.01) +
    scale_x_discrete(labels = treatment_ionis_rename,
                     guide = guide_axis(angle = 45)) +
    scale_y_continuous(breaks = pretty_breaks()) +
    scale_colour_manual(values = treatment_ionis_colour) +
    scale_fill_manual(values = treatment_ionis_colour) +
    labs(x = "Treatment",
         y = "Relative expression (%)",
         colour = "Treatment",
         fill = "Treatment") +
    theme_minimal() +
    theme(text = element_text(size = element_text_size),
          legend.position = "none")
  
  # Output
  out <- mget(c("ttest", "my_pwc", "my_plot"))
  return(out)
  
},
plot_data = plot_data,
padj_method = padj_method,
hide_ns = hide_ns,
treatment_ionis_rename = treatment_ionis_rename,
treatment_ionis_colour = treatment_ionis_colour,
element_text_size = element_text_size)

# Display ttests
my_ttests <- lapply(my_experiment_names, function(my_experiment_name,
                                                  my_plots) {
  my_plots[[my_experiment_name]]$ttest
},
my_plots = my_plots)
for (my_experiment_name in my_experiment_names) {
  my_ttest <- my_ttests[[my_experiment_name]]
  print(kable(my_ttest, caption = paste(my_experiment_name, "t-test")))
}
## 
## 
## Table: invivo_titration t-test
## 
## |tissue      |.y.   |group1  |group2 | n1| n2|  statistic|       df|        p|     p.adj|p.adj.signif | y.position|groups        |
## |:-----------|:-----|:-------|:------|--:|--:|----------:|--------:|--------:|---------:|:------------|----------:|:-------------|
## |Cortex      |level |Vehicle |3ug    |  4|  4|  1.8542000| 3.934168| 1.39e-01| 0.1516364|ns           |   116.4359|Vehicle, 3ug  |
## |Cortex      |level |Vehicle |10ug   |  4|  4|  1.9604778| 3.502115| 1.32e-01| 0.1508571|ns           |   128.4561|Vehicle, 10ug |
## |Cortex      |level |Vehicle |30ug   |  4|  5| 12.0544775| 6.802275| 7.70e-06| 0.0000449|***          |   140.4764|Vehicle, 30ug |
## |Cortex      |level |3ug     |10ug   |  4|  4|  0.4778913| 5.476986| 6.51e-01| 0.6510000|ns           |   152.4966|3ug , 10ug    |
## |Cortex      |level |3ug     |30ug   |  4|  5|  5.4406421| 4.771674| 3.00e-03| 0.0055385|**           |   164.5169|3ug , 30ug    |
## |Cortex      |level |10ug    |30ug   |  4|  5|  3.6419864| 3.964163| 2.20e-02| 0.0368000|*            |   176.5372|10ug, 30ug    |
## |Striatum    |level |Vehicle |3ug    |  4|  4|  0.5318878| 4.900537| 6.18e-01| 0.6448696|ns           |   120.4869|Vehicle, 3ug  |
## |Striatum    |level |Vehicle |10ug   |  4|  4|  2.2498355| 3.606236| 9.50e-02| 0.1200000|ns           |   132.5071|Vehicle, 10ug |
## |Striatum    |level |Vehicle |30ug   |  4|  5|  7.1738608| 4.011652| 2.00e-03| 0.0040000|**           |   144.5274|Vehicle, 30ug |
## |Striatum    |level |3ug     |10ug   |  4|  4|  2.5708958| 4.585510| 5.40e-02| 0.0720000|ns           |   156.5476|3ug , 10ug    |
## |Striatum    |level |3ug     |30ug   |  4|  5|  9.8404489| 5.570527| 9.93e-05| 0.0002648|***          |   168.5679|3ug , 30ug    |
## |Striatum    |level |10ug    |30ug   |  4|  5| 10.3652603| 6.925037| 1.82e-05| 0.0000624|***          |   180.5882|10ug, 30ug    |
## |Brainstem   |level |Vehicle |3ug    |  4|  4|  6.2762166| 5.062854| 1.00e-03| 0.0021818|**           |   113.6829|Vehicle, 3ug  |
## |Brainstem   |level |Vehicle |10ug   |  4|  4| 17.6634697| 5.767722| 3.00e-06| 0.0000244|***          |   125.7031|Vehicle, 10ug |
## |Brainstem   |level |Vehicle |30ug   |  4|  5| 32.4960641| 6.975689| 0.00e+00| 0.0000002|***          |   137.7234|Vehicle, 30ug |
## |Brainstem   |level |3ug     |10ug   |  4|  4|  8.0782464| 5.643734| 2.62e-04| 0.0006288|***          |   149.7436|3ug , 10ug    |
## |Brainstem   |level |3ug     |30ug   |  4|  5| 18.8693737| 5.569056| 2.90e-06| 0.0000244|***          |   161.7639|3ug , 30ug    |
## |Brainstem   |level |10ug    |30ug   |  4|  5| 12.2559102| 6.544504| 9.40e-06| 0.0000449|***          |   173.7842|10ug, 30ug    |
## |Spinal cord |level |Vehicle |3ug    |  4|  4|  1.8137090| 5.196380| 1.27e-01| 0.1508571|ns           |   114.0809|Vehicle, 3ug  |
## |Spinal cord |level |Vehicle |10ug   |  4|  4|  4.1886198| 3.144929| 2.30e-02| 0.0368000|*            |   126.1011|Vehicle, 10ug |
## |Spinal cord |level |Vehicle |30ug   |  4|  5| 15.5784053| 4.806628| 2.68e-05| 0.0000804|***          |   138.1214|Vehicle, 30ug |
## |Spinal cord |level |3ug     |10ug   |  4|  4|  3.6277004| 3.331969| 3.00e-02| 0.0423529|*            |   150.1416|3ug , 10ug    |
## |Spinal cord |level |3ug     |30ug   |  4|  5| 13.7512712| 5.706043| 1.37e-05| 0.0000548|***          |   162.1619|3ug , 30ug    |
## |Spinal cord |level |10ug    |30ug   |  4|  5|  3.2870097| 4.377864| 2.60e-02| 0.0390000|*            |   174.1822|10ug, 30ug    |
## 
## 
## Table: invivo_titration_SCR t-test
## 
## |tissue      |.y.   |group1  |group2 | n1| n2|  statistic|       df|     p|     p.adj|p.adj.signif | y.position|groups          |
## |:-----------|:-----|:-------|:------|--:|--:|----------:|--------:|-----:|---------:|:------------|----------:|:---------------|
## |Cortex      |level |Vehicle |3ug    |  4|  4|  0.1049662| 5.594147| 0.920| 0.9200000|ns           |   115.1301|Vehicle, 3ug    |
## |Cortex      |level |Vehicle |10ug   |  4|  4|  2.2245897| 5.918406| 0.068| 0.1133333|ns           |   121.5746|Vehicle, 10ug   |
## |Cortex      |level |Vehicle |30ug   |  4|  3|  6.7692358| 2.581470| 0.010| 0.0500000|*            |   128.0191|Vehicle, 30ug   |
## |Cortex      |level |Vehicle |SCRaso |  4|  4| -0.8455827| 5.916147| 0.431| 0.4600000|ns           |   134.4637|Vehicle, SCRaso |
## |Cortex      |level |3ug     |10ug   |  4|  4|  1.8321285| 5.855849| 0.118| 0.1653333|ns           |   140.9082|3ug , 10ug      |
## |Cortex      |level |3ug     |30ug   |  4|  3|  6.4166254| 3.000136| 0.008| 0.0500000|*            |   147.3527|3ug , 30ug      |
## |Cortex      |level |3ug     |SCRaso |  4|  4| -0.8348014| 5.858741| 0.437| 0.4600000|ns           |   153.7973|3ug   , SCRaso  |
## |Cortex      |level |10ug    |30ug   |  4|  3|  5.4902697| 2.734738| 0.015| 0.0500000|*            |   160.2418|10ug, 30ug      |
## |Cortex      |level |10ug    |SCRaso |  4|  4| -2.9027231| 5.999983| 0.027| 0.0700000|ns           |   166.6863|10ug  , SCRaso  |
## |Cortex      |level |30ug    |SCRaso |  3|  4| -7.1014215| 2.737148| 0.008| 0.0500000|*            |   173.1309|30ug  , SCRaso  |
## |Spinal cord |level |Vehicle |3ug    |  4|  4|  3.5857849| 5.441242| 0.014| 0.0500000|*            |   113.0961|Vehicle, 3ug    |
## |Spinal cord |level |Vehicle |10ug   |  4|  4|  4.4059859| 4.737829| 0.008| 0.0500000|*            |   119.5406|Vehicle, 10ug   |
## |Spinal cord |level |Vehicle |30ug   |  4|  3|  4.8340610| 2.218758| 0.033| 0.0700000|ns           |   125.9851|Vehicle, 30ug   |
## |Spinal cord |level |Vehicle |SCRaso |  4|  4|  1.2633565| 5.058685| 0.262| 0.3082353|ns           |   132.4297|Vehicle, SCRaso |
## |Spinal cord |level |3ug     |10ug   |  4|  4|  1.2454440| 5.687328| 0.262| 0.3082353|ns           |   138.8742|3ug , 10ug      |
## |Spinal cord |level |3ug     |30ug   |  4|  3|  3.3573428| 2.426350| 0.060| 0.1090909|ns           |   145.3187|3ug , 30ug      |
## |Spinal cord |level |3ug     |SCRaso |  4|  4| -1.7908931| 5.902224| 0.124| 0.1653333|ns           |   151.7633|3ug   , SCRaso  |
## |Spinal cord |level |10ug    |30ug   |  4|  3|  2.6587774| 2.686702| 0.086| 0.1323077|ns           |   158.2078|10ug, 30ug      |
## |Spinal cord |level |10ug    |SCRaso |  4|  4| -2.7716089| 5.929482| 0.033| 0.0700000|ns           |   164.6523|10ug  , SCRaso  |
## |Spinal cord |level |30ug    |SCRaso |  3|  4| -4.1376927| 2.552372| 0.035| 0.0700000|ns           |   171.0969|30ug  , SCRaso  |
# Display ANOVA posthoc
my_pwcs <- lapply(my_experiment_names, function(my_experiment_name,
                                                my_plots) {
  my_plots[[my_experiment_name]]$my_pwc
},
my_plots = my_plots)
for (my_experiment_name in my_experiment_names) {
  my_pwc <- my_pwcs[[my_experiment_name]]
  print(kable(my_pwc, caption = paste(my_experiment_name, "ANOVA with corrected posthoc pairwise comparisons")))
}
## 
## 
## Table: invivo_titration ANOVA with corrected posthoc pairwise comparisons
## 
## |contrast       |group1  |group2 |tissue      |  estimate|       SE| df|    t.ratio|   p.value|     p.adj|p.adj.signif |
## |:--------------|:-------|:------|:-----------|---------:|--------:|--:|----------:|---------:|---------:|:------------|
## |Vehicle - 3ug  |Vehicle |3ug    |Cortex      | 12.985657| 6.591570| 52|  1.9700402| 0.2125051| 0.2550061|ns           |
## |Vehicle - 10ug |Vehicle |10ug   |Cortex      | 18.273183| 6.591570| 52|  2.7722052| 0.0374764| 0.0499685|*            |
## |Vehicle - 30ug |Vehicle |30ug   |Cortex      | 53.393389| 6.253312| 52|  8.5384178| 0.0000000| 0.0000000|***          |
## |3ug - 10ug     |3ug     |10ug   |Cortex      |  5.287527| 6.591570| 52|  0.8021650| 0.8530927| 0.9016176|ns           |
## |3ug - 30ug     |3ug     |30ug   |Cortex      | 40.407732| 6.253312| 52|  6.4618131| 0.0000002| 0.0000005|***          |
## |10ug - 30ug    |10ug    |30ug   |Cortex      | 35.120206| 6.253312| 52|  5.6162569| 0.0000045| 0.0000084|***          |
## |Vehicle - 3ug  |Vehicle |3ug    |Striatum    |  3.659846| 6.591570| 52|  0.5552313| 0.9447105| 0.9447105|ns           |
## |Vehicle - 10ug |Vehicle |10ug   |Striatum    | 13.950394| 6.591570| 52|  2.1163994| 0.1614185| 0.2038971|ns           |
## |Vehicle - 30ug |Vehicle |30ug   |Striatum    | 45.805947| 6.253312| 52|  7.3250699| 0.0000000| 0.0000000|***          |
## |3ug - 10ug     |3ug     |10ug   |Striatum    | 10.290548| 6.591570| 52|  1.5611681| 0.4094088| 0.4678958|ns           |
## |3ug - 30ug     |3ug     |30ug   |Striatum    | 42.146101| 6.253312| 52|  6.7398047| 0.0000001| 0.0000002|***          |
## |10ug - 30ug    |10ug    |30ug   |Striatum    | 31.855553| 6.253312| 52|  5.0941891| 0.0000288| 0.0000461|***          |
## |Vehicle - 3ug  |Vehicle |3ug    |Brainstem   | 18.896578| 6.591570| 52|  2.8667797| 0.0295351| 0.0416967|*            |
## |Vehicle - 10ug |Vehicle |10ug   |Brainstem   | 44.900223| 6.591570| 52|  6.8117651| 0.0000001| 0.0000002|***          |
## |Vehicle - 30ug |Vehicle |30ug   |Brainstem   | 77.182277| 6.253312| 52| 12.3426240| 0.0000000| 0.0000000|***          |
## |3ug - 10ug     |3ug     |10ug   |Brainstem   | 26.003645| 6.591570| 52|  3.9449854| 0.0013358| 0.0020037|**           |
## |3ug - 30ug     |3ug     |30ug   |Brainstem   | 58.285700| 6.253312| 52|  9.3207729| 0.0000000| 0.0000000|***          |
## |10ug - 30ug    |10ug    |30ug   |Brainstem   | 32.282054| 6.253312| 52|  5.1623931| 0.0000227| 0.0000389|***          |
## |Vehicle - 3ug  |Vehicle |3ug    |Spinal cord |  5.128342| 6.591570| 52|  0.7780154| 0.8640502| 0.9016176|ns           |
## |Vehicle - 10ug |Vehicle |10ug   |Spinal cord | 42.464438| 6.591570| 52|  6.4422348| 0.0000002| 0.0000005|***          |
## |Vehicle - 30ug |Vehicle |30ug   |Spinal cord | 79.016220| 6.253312| 52| 12.6358995| 0.0000000| 0.0000000|***          |
## |3ug - 10ug     |3ug     |10ug   |Spinal cord | 37.336096| 6.591570| 52|  5.6642194| 0.0000038| 0.0000076|***          |
## |3ug - 30ug     |3ug     |30ug   |Spinal cord | 73.887878| 6.253312| 52| 11.8157992| 0.0000000| 0.0000000|***          |
## |10ug - 30ug    |10ug    |30ug   |Spinal cord | 36.551782| 6.253312| 52|  5.8451877| 0.0000020| 0.0000043|***          |
## 
## 
## Table: invivo_titration_SCR ANOVA with corrected posthoc pairwise comparisons
## 
## |contrast         |group1  |group2 |tissue      |    estimate|       SE| df|    t.ratio|   p.value|     p.adj|p.adj.signif |
## |:----------------|:-------|:------|:-----------|-----------:|--------:|--:|----------:|---------:|---------:|:------------|
## |Vehicle - 3ug    |Vehicle |3ug    |Cortex      |   0.3914004| 6.702412| 28|  0.0583970| 0.9999972| 0.9999972|ns           |
## |Vehicle - 10ug   |Vehicle |10ug   |Cortex      |   7.5474167| 6.702412| 28|  1.1260747| 0.7915231| 0.9661809|ns           |
## |Vehicle - 30ug   |Vehicle |30ug   |Cortex      |  43.0115304| 7.239432| 28|  5.9412850| 0.0000200| 0.0000771|***          |
## |Vehicle - SCRaso |Vehicle |SCRaso |Cortex      |  -2.8714926| 6.702412| 28| -0.4284268| 0.9925764| 0.9999972|ns           |
## |3ug - 10ug       |3ug     |10ug   |Cortex      |   7.1560163| 6.702412| 28|  1.0676778| 0.8212538| 0.9661809|ns           |
## |3ug - 30ug       |3ug     |30ug   |Cortex      |  42.6201300| 7.239432| 28|  5.8872199| 0.0000231| 0.0000771|***          |
## |3ug - SCRaso     |3ug     |SCRaso |Cortex      |  -3.2628930| 6.702412| 28| -0.4868237| 0.9879628| 0.9999972|ns           |
## |10ug - 30ug      |10ug    |30ug   |Cortex      |  35.4641137| 7.239432| 28|  4.8987424| 0.0003293| 0.0008232|***          |
## |10ug - SCRaso    |10ug    |SCRaso |Cortex      | -10.4189093| 6.702412| 28| -1.5545015| 0.5374236| 0.8268055|ns           |
## |30ug - SCRaso    |30ug    |SCRaso |Cortex      | -45.8830230| 7.239432| 28| -6.3379311| 0.0000070| 0.0000349|***          |
## |Vehicle - 3ug    |Vehicle |3ug    |Spinal cord |  18.7403032| 6.702412| 28|  2.7960536| 0.0646660| 0.1175745|ns           |
## |Vehicle - 10ug   |Vehicle |10ug   |Spinal cord |  27.2889479| 6.702412| 28|  4.0715116| 0.0029549| 0.0065663|**           |
## |Vehicle - 30ug   |Vehicle |30ug   |Spinal cord |  64.8874611| 7.239432| 28|  8.9630593| 0.0000000| 0.0000002|***          |
## |Vehicle - SCRaso |Vehicle |SCRaso |Spinal cord |   7.2179372| 6.702412| 28|  1.0769164| 0.8166817| 0.9661809|ns           |
## |3ug - 10ug       |3ug     |10ug   |Spinal cord |   8.5486447| 6.702412| 28|  1.2754580| 0.7078962| 0.9661809|ns           |
## |3ug - 30ug       |3ug     |30ug   |Spinal cord |  46.1471579| 7.239432| 28|  6.3744167| 0.0000063| 0.0000349|***          |
## |3ug - SCRaso     |3ug     |SCRaso |Spinal cord | -11.5223660| 6.702412| 28| -1.7191372| 0.4390542| 0.7317570|ns           |
## |10ug - 30ug      |10ug    |30ug   |Spinal cord |  37.5985132| 7.239432| 28|  5.1935720| 0.0001490| 0.0004257|***          |
## |10ug - SCRaso    |10ug    |SCRaso |Spinal cord | -20.0710107| 6.702412| 28| -2.9945952| 0.0416897| 0.0833794|ns           |
## |30ug - SCRaso    |30ug    |SCRaso |Spinal cord | -57.6695239| 7.239432| 28| -7.9660285| 0.0000001| 0.0000011|***          |
# Display plots
my_bars <- lapply(my_experiment_names, function(my_experiment_name,
                                                my_plots) {
  my_plots[[my_experiment_name]]$my_plot
},
my_plots = my_plots)
for (my_experiment_name in my_experiment_names) {
  print(my_experiment_name)
  print(my_bars[[my_experiment_name]])
}
## [1] "invivo_titration"

## [1] "invivo_titration_SCR"

# Export plots
# for (my_experiment_name in my_experiment_names) {
#   for (output_format in output_formats) {
#   ggsave(filename = file.path(out_dir, "in_vivo_titration", paste0("in_vivo_titration", "_", my_experiment_name, ".", output_format)),
#          plot = my_bars[[my_experiment_name]],
#          height = 6, width = 12)
#   }
# }

# Export experiment 1 plot
for (output_format in output_formats) {
  ggsave(filename = file.path(out_dir, "in_vivo_titration", paste0("in_vivo_titration-experiment1", ".", output_format)),
         plot = my_bars[[my_experiment_names[[1]]]],
         height = 6, width = 12)
}

# Export experiment 2 plot
for (output_format in output_formats) {
  ggsave(filename = file.path(out_dir, "in_vivo_titration", paste0("in_vivo_titration-experiment2", ".", output_format)),
         plot = my_bars[[my_experiment_names[[2]]]],
         height = 6, width = 9)
}

Combined experiments with SCR separate

# List of plots
my_treatment_groups <- list(MSH3aso = c("Vehicle", "3ug", "10ug", "30ug"),
                            SCRaso = c("Vehicle", "SCRaso"))

# Plot
my_plots <- lapply(my_treatment_groups, function(my_treatment_group,
                                                 plot_data,
                                                 padj_method,
                                                 hide_ns,
                                                 treatment_ionis_rename,
                                                 treatment_ionis_colour,
                                                 element_text_size) {
  
  # Filter data
  my_data <- plot_data %>%
    dplyr::filter(treatment %in% my_treatment_group)
  
  # Filter to keep only tissues where all treatments are present
  my_data <- my_data %>%
    group_by(tissue) %>%
    filter(all(my_treatment_group %in% treatment)) %>%
    ungroup()
  
  # t-test
  ttest <- 
    my_data %>%
    group_by(tissue) %>%
    t_test(level ~ treatment) %>%
    adjust_pvalue(method = padj_method) %>%
    add_significance(cutpoints = c(0, 0.001, 0.01, 0.05, 1),
                     symbols = c("***", "**", "*", "ns")) %>%
    add_y_position() %>%
    ungroup()
  
  
  # ANOVA and corrected posthoc pairwise comparisons
  my_aov <- aov(level ~ treatment * tissue, data = my_data)
  my_emmeans <- emmeans(my_aov, pairwise ~ treatment | tissue) # emmeans 'adjust' argument is giving odd results, with some p values getting more significant
  my_pwc <- as.data.frame(my_emmeans$contrasts) %>%
    dplyr::mutate(group1 = str_split_fixed(contrast, " - ", 2)[, 1],
                  group2 = str_split_fixed(contrast, " - ", 2)[, 2],
                  p.adj = p.adjust(p.value, method = padj_method),
                  p.adj.signif = symnum(p.adj, corr = FALSE, na = FALSE, legend = FALSE,
                                        cutpoints = c(0, 0.001, 0.01, 0.05, 1), 
                                        symbols = c("***", "**", "*", "ns"))) %>%
    relocate(c(group1, group2), .after = "contrast") %>%
    relocate(c(p.adj, p.adj.signif), .after = "p.value")
  if (hide_ns) {
    my_pwc_plot <- my_pwc %>%
      dplyr::filter(p.adj < 0.05)
  } else {
    my_pwc_plot <- my_pwc
  }
  y_max <- my_data %>%
    group_by(tissue) %>%
    dplyr::summarise(ymax = max(level, na.rm = TRUE))
  if (nrow(my_pwc_plot) > 0) {
    my_pwc_plot <- my_pwc_plot %>%
      dplyr::left_join(y_max, by = join_by(tissue)) %>%
      group_by(tissue) %>%
      dplyr::mutate(y.position = ymax + seq(0.1 * unique(ymax), 0.5 * unique(ymax), length.out = n())) %>%
      ungroup()
  } else {
    my_pwc_plot[1,] <- NA
    my_pwc_plot$y.position <- NA
  }
  
  # Plot
  my_plot <-
    ggplot(my_data, aes(x = treatment, y = level)) +
    facet_wrap(vars(tissue), nrow = 1, scales = "free_x") +
    stat_summary(aes(fill = treatment),
                 # fun = mean, colour = "black", geom = "bar", alpha = 0.8) +
                 fun = mean, colour = "black", geom = "bar") +
    stat_summary(fun.data = mean_se,  
                 geom = "errorbar", width = 0.2) +
    geom_jitter(aes(fill = treatment),
                width = 0.2, height = 0, size = point_size, colour = "black", shape = 21) +
    geom_hline(yintercept = 0, colour = "darkgrey") +
    geom_hline(yintercept = 100, colour = "darkgrey", linetype = "dashed") +
    # stat_pvalue_manual(ttest, hide.ns = hide_ns, label = "p.adj.signif")
    stat_pvalue_manual(my_pwc_plot, label = "p.adj.signif", hide.ns = hide_ns, tip.length = 0.01) +
    scale_x_discrete(labels = treatment_ionis_rename,
                     guide = guide_axis(angle = 45)) +
    scale_y_continuous(breaks = pretty_breaks()) +
    scale_colour_manual(values = treatment_ionis_colour) +
    scale_fill_manual(values = treatment_ionis_colour) +
    labs(x = "Treatment",
         y = "Relative expression (%)",
         colour = "Treatment",
         fill = "Treatment") +
    theme_minimal() +
    theme(text = element_text(size = element_text_size),
          legend.position = "none")
  
  # Output
  out <- mget(c("ttest", "my_pwc", "my_plot"))
  return(out)
  
},
plot_data = plot_data,
padj_method = padj_method,
hide_ns = hide_ns,
treatment_ionis_rename = treatment_ionis_rename,
treatment_ionis_colour = treatment_ionis_colour,
element_text_size = element_text_size)


# Vector of treatment group names
my_treatment_group_names <- setNames(names(my_treatment_groups), names(my_treatment_groups))


# Display ttests
my_ttests <- lapply(my_treatment_group_names, function(my_treatment_group_name,
                                                       my_plots) {
  my_plots[[my_treatment_group_name]]$ttest
},
my_plots = my_plots)
for (my_treatment_group_name in my_treatment_group_names) {
  my_ttest <- my_ttests[[my_treatment_group_name]]
  print(kable(my_ttest, caption = paste(my_treatment_group_name, "t-test")))
}
## 
## 
## Table: MSH3aso t-test
## 
## |tissue      |.y.   |group1  |group2 | n1| n2|  statistic|        df|        p|     p.adj|p.adj.signif | y.position|groups        |
## |:-----------|:-----|:-------|:------|--:|--:|----------:|---------:|--------:|---------:|:------------|----------:|:-------------|
## |Cortex      |level |Vehicle |3ug    |  8|  8|  1.5415216|  9.132123| 1.57e-01| 0.1712727|ns           |   116.4359|Vehicle, 3ug  |
## |Cortex      |level |Vehicle |10ug   |  8|  8|  2.5865314|  8.569527| 3.00e-02| 0.0378947|*            |   128.4561|Vehicle, 10ug |
## |Cortex      |level |Vehicle |30ug   |  8|  8| 13.0916812|  9.896674| 1.00e-07| 0.0000017|***          |   140.4764|Vehicle, 30ug |
## |Cortex      |level |3ug     |10ug   |  8|  8|  1.0006693| 13.661257| 3.34e-01| 0.3485217|ns           |   152.4966|3ug , 10ug    |
## |Cortex      |level |3ug     |30ug   |  8|  8|  8.0797152| 13.639093| 1.50e-06| 0.0000079|***          |   164.5169|3ug , 30ug    |
## |Cortex      |level |10ug    |30ug   |  8|  8|  6.2577056| 12.757057| 3.20e-05| 0.0000768|***          |   176.5372|10ug, 30ug    |
## |Striatum    |level |Vehicle |3ug    |  4|  4|  0.5318878|  4.900537| 6.18e-01| 0.6180000|ns           |   120.4869|Vehicle, 3ug  |
## |Striatum    |level |Vehicle |10ug   |  4|  4|  2.2498355|  3.606236| 9.50e-02| 0.1085714|ns           |   132.5071|Vehicle, 10ug |
## |Striatum    |level |Vehicle |30ug   |  4|  5|  7.1738608|  4.011652| 2.00e-03| 0.0030000|**           |   144.5274|Vehicle, 30ug |
## |Striatum    |level |3ug     |10ug   |  4|  4|  2.5708958|  4.585510| 5.40e-02| 0.0648000|ns           |   156.5476|3ug , 10ug    |
## |Striatum    |level |3ug     |30ug   |  4|  5|  9.8404489|  5.570527| 9.93e-05| 0.0002167|***          |   168.5679|3ug , 30ug    |
## |Striatum    |level |10ug    |30ug   |  4|  5| 10.3652603|  6.925037| 1.82e-05| 0.0000485|***          |   180.5882|10ug, 30ug    |
## |Brainstem   |level |Vehicle |3ug    |  4|  4|  6.2762166|  5.062854| 1.00e-03| 0.0016000|**           |   113.6829|Vehicle, 3ug  |
## |Brainstem   |level |Vehicle |10ug   |  4|  4| 17.6634697|  5.767722| 3.00e-06| 0.0000105|***          |   125.7031|Vehicle, 10ug |
## |Brainstem   |level |Vehicle |30ug   |  4|  5| 32.4960641|  6.975689| 0.00e+00| 0.0000002|***          |   137.7234|Vehicle, 30ug |
## |Brainstem   |level |3ug     |10ug   |  4|  4|  8.0782464|  5.643734| 2.62e-04| 0.0005240|***          |   149.7436|3ug , 10ug    |
## |Brainstem   |level |3ug     |30ug   |  4|  5| 18.8693737|  5.569056| 2.90e-06| 0.0000105|***          |   161.7639|3ug , 30ug    |
## |Brainstem   |level |10ug    |30ug   |  4|  5| 12.2559102|  6.544504| 9.40e-06| 0.0000281|***          |   173.7842|10ug, 30ug    |
## |Spinal cord |level |Vehicle |3ug    |  8|  8|  3.1598645|  9.857447| 1.00e-02| 0.0133333|*            |   117.3129|Vehicle, 3ug  |
## |Spinal cord |level |Vehicle |10ug   |  8|  8|  5.6296674|  7.974333| 4.99e-04| 0.0008554|***          |   129.3331|Vehicle, 10ug |
## |Spinal cord |level |Vehicle |30ug   |  8|  8| 12.3280829|  8.050690| 1.70e-06| 0.0000079|***          |   141.3534|Vehicle, 30ug |
## |Spinal cord |level |3ug     |10ug   |  8|  8|  3.3243802| 11.143146| 7.00e-03| 0.0098824|**           |   153.3736|3ug , 10ug    |
## |Spinal cord |level |3ug     |30ug   |  8|  8|  9.2092605| 11.401101| 1.30e-06| 0.0000079|***          |   165.3939|3ug , 30ug    |
## |Spinal cord |level |10ug    |30ug   |  8|  8|  4.6709975| 13.979700| 3.62e-04| 0.0006683|***          |   177.4142|10ug, 30ug    |
## 
## 
## Table: SCRaso t-test
## 
## |tissue      |.y.   |group1  |group2 | n1| n2|  statistic|       df|     p| p.adj|p.adj.signif | y.position|groups          |
## |:-----------|:-----|:-------|:------|--:|--:|----------:|--------:|-----:|-----:|:------------|----------:|:---------------|
## |Cortex      |level |Vehicle |SCRaso |  8|  4| -0.8948261| 5.464449| 0.409| 0.409|ns           |   110.1899|Vehicle, SCRaso |
## |Spinal cord |level |Vehicle |SCRaso |  8|  4|  1.4272505| 3.660822| 0.233| 0.409|ns           |   108.1559|Vehicle, SCRaso |
# Display ANOVA posthoc
my_pwcs <- lapply(my_treatment_group_names, function(my_treatment_group_name,
                                                     my_plots) {
  my_plots[[my_treatment_group_name]]$my_pwc
},
my_plots = my_plots)
for (my_treatment_group_name in my_treatment_group_names) {
  my_pwc <- my_pwcs[[my_treatment_group_name]]
  print(kable(my_pwc, caption = paste(my_treatment_group_name, "ANOVA with corrected posthoc pairwise comparisons")))
}
## 
## 
## Table: MSH3aso ANOVA with corrected posthoc pairwise comparisons
## 
## |contrast       |group1  |group2 |tissue      |  estimate|       SE| df|    t.ratio|   p.value|     p.adj|p.adj.signif |
## |:--------------|:-------|:------|:-----------|---------:|--------:|--:|----------:|---------:|---------:|:------------|
## |Vehicle - 3ug  |Vehicle |3ug    |Cortex      |  6.688529| 5.113932| 82|  1.3079033| 0.5605969| 0.6115603|ns           |
## |Vehicle - 10ug |Vehicle |10ug   |Cortex      | 12.910300| 5.113932| 82|  2.5245348| 0.0635297| 0.0847063|ns           |
## |Vehicle - 30ug |Vehicle |30ug   |Cortex      | 49.453178| 5.113932| 82|  9.6702841| 0.0000000| 0.0000000|***          |
## |3ug - 10ug     |3ug     |10ug   |Cortex      |  6.221771| 5.113932| 82|  1.2166316| 0.6181520| 0.6450282|ns           |
## |3ug - 30ug     |3ug     |30ug   |Cortex      | 42.764649| 5.113932| 82|  8.3623808| 0.0000000| 0.0000000|***          |
## |10ug - 30ug    |10ug    |30ug   |Cortex      | 36.542878| 5.113932| 82|  7.1457493| 0.0000000| 0.0000000|***          |
## |Vehicle - 3ug  |Vehicle |3ug    |Striatum    |  3.659846| 7.232192| 82|  0.5060493| 0.9574291| 0.9574291|ns           |
## |Vehicle - 10ug |Vehicle |10ug   |Striatum    | 13.950394| 7.232192| 82|  1.9289301| 0.2240923| 0.2689108|ns           |
## |Vehicle - 30ug |Vehicle |30ug   |Striatum    | 45.805947| 6.861060| 82|  6.6762200| 0.0000000| 0.0000000|***          |
## |3ug - 10ug     |3ug     |10ug   |Striatum    | 10.290548| 7.232192| 82|  1.4228808| 0.4888740| 0.5587132|ns           |
## |3ug - 30ug     |3ug     |30ug   |Striatum    | 42.146101| 6.861060| 82|  6.1427972| 0.0000002| 0.0000003|***          |
## |10ug - 30ug    |10ug    |30ug   |Striatum    | 31.855553| 6.861060| 82|  4.6429491| 0.0000752| 0.0001289|***          |
## |Vehicle - 3ug  |Vehicle |3ug    |Brainstem   | 18.896578| 7.232192| 82|  2.6128422| 0.0512146| 0.0723029|ns           |
## |Vehicle - 10ug |Vehicle |10ug   |Brainstem   | 44.900223| 7.232192| 82|  6.2083834| 0.0000001| 0.0000003|***          |
## |Vehicle - 30ug |Vehicle |30ug   |Brainstem   | 77.182277| 6.861060| 82| 11.2493224| 0.0000000| 0.0000000|***          |
## |3ug - 10ug     |3ug     |10ug   |Brainstem   | 26.003645| 7.232192| 82|  3.5955412| 0.0030423| 0.0045635|**           |
## |3ug - 30ug     |3ug     |30ug   |Brainstem   | 58.285700| 6.861060| 82|  8.4951449| 0.0000000| 0.0000000|***          |
## |10ug - 30ug    |10ug    |30ug   |Brainstem   | 32.282054| 6.861060| 82|  4.7051117| 0.0000593| 0.0001095|***          |
## |Vehicle - 3ug  |Vehicle |3ug    |Spinal cord | 11.934323| 5.113932| 82|  2.3336881| 0.0987803| 0.1247751|ns           |
## |Vehicle - 10ug |Vehicle |10ug   |Spinal cord | 34.876693| 5.113932| 82|  6.8199365| 0.0000000| 0.0000000|***          |
## |Vehicle - 30ug |Vehicle |30ug   |Spinal cord | 73.707582| 5.113932| 82| 14.4130931| 0.0000000| 0.0000000|***          |
## |3ug - 10ug     |3ug     |10ug   |Spinal cord | 22.942370| 5.113932| 82|  4.4862484| 0.0001357| 0.0002171|***          |
## |3ug - 30ug     |3ug     |30ug   |Spinal cord | 61.773259| 5.113932| 82| 12.0794051| 0.0000000| 0.0000000|***          |
## |10ug - 30ug    |10ug    |30ug   |Spinal cord | 38.830889| 5.113932| 82|  7.5931567| 0.0000000| 0.0000000|***          |
## 
## 
## Table: SCRaso ANOVA with corrected posthoc pairwise comparisons
## 
## |contrast         |group1  |group2 |tissue      |  estimate|       SE| df|    t.ratio|   p.value|     p.adj|p.adj.signif |
## |:----------------|:-------|:------|:-----------|---------:|--------:|--:|----------:|---------:|---------:|:------------|
## |Vehicle - SCRaso |Vehicle |SCRaso |Cortex      | -2.683435| 3.465514| 20| -0.7743251| 0.4477962| 0.4477962|ns           |
## |Vehicle - SCRaso |Vehicle |SCRaso |Spinal cord |  7.259351| 3.465514| 20|  2.0947398| 0.0491315| 0.0982631|ns           |
# Display plots
my_bars <- lapply(my_treatment_group_names, function(my_treatment_group_name,
                                                     my_plots) {
  my_plots[[my_treatment_group_name]]$my_plot
},
my_plots = my_plots)
for (my_treatment_group_name in my_treatment_group_names) {
  print(my_treatment_group_name)
  print(my_bars[[my_treatment_group_name]])
}
## [1] "MSH3aso"

## [1] "SCRaso"

# Export plots
# for (my_experiment_name in my_experiment_names) {
#   for (output_format in output_formats) {
#   ggsave(filename = file.path(out_dir, "in_vivo_titration", paste0("in_vivo_titration", "_", my_experiment_name, ".", output_format)),
#          plot = my_bars[[my_experiment_name]],
#          height = 6, width = 12)
#   }
# }

# Export experiment 1 plot
for (output_format in output_formats) {
  ggsave(filename = file.path(out_dir, "in_vivo_titration", paste0("in_vivo_titration-combined-", my_treatment_group_names[[1]], ".", output_format)),
         plot = my_bars[[my_treatment_group_names[[1]]]],
         height = 6, width = 12)
}

# Export experiment 2 plot
for (output_format in output_formats) {
  ggsave(filename = file.path(out_dir, "in_vivo_titration", paste0("in_vivo_titration-combined-", my_treatment_group_names[[2]], ".", output_format)),
         plot = my_bars[[my_treatment_group_names[[2]]]],
         height = 5, width = 7)
}

Custom

Experiment 2 with only Vehicle, SCR ASO and the highest dose of MSH3 ASO.

# Filter data
my_data <- plot_data %>%
  dplyr::filter(experiment_id == "invivo_titration_SCR",
                treatment %in% c("Vehicle", "SCRaso", "30ug"))

# t-test
ttest <- 
  my_data %>%
  group_by(tissue) %>%
  t_test(level ~ treatment) %>%
  adjust_pvalue(method = padj_method) %>%
  add_significance(cutpoints = c(0, 0.001, 0.01, 0.05, 1),
                   symbols = c("***", "**", "*", "ns")) %>%
  add_y_position() %>%
  ungroup()

# Display ttest
print(kable(ttest, caption = "t-test"))
## 
## 
## Table: t-test
## 
## |tissue      |.y.   |group1  |group2 | n1| n2|  statistic|       df|     p|  p.adj|p.adj.signif | y.position|groups          |
## |:-----------|:-----|:-------|:------|--:|--:|----------:|--------:|-----:|------:|:------------|----------:|:---------------|
## |Cortex      |level |Vehicle |30ug   |  4|  3|  6.7692358| 2.581470| 0.010| 0.0300|*            |   115.1301|Vehicle, 30ug   |
## |Cortex      |level |Vehicle |SCRaso |  4|  4| -0.8455827| 5.916147| 0.431| 0.4310|ns           |   123.8302|Vehicle, SCRaso |
## |Cortex      |level |30ug    |SCRaso |  3|  4| -7.1014215| 2.737148| 0.008| 0.0300|*            |   132.5303|30ug  , SCRaso  |
## |Spinal cord |level |Vehicle |30ug   |  4|  3|  4.8340610| 2.218758| 0.033| 0.0525|ns           |   113.0961|Vehicle, 30ug   |
## |Spinal cord |level |Vehicle |SCRaso |  4|  4|  1.2633565| 5.058685| 0.262| 0.3144|ns           |   121.7962|Vehicle, SCRaso |
## |Spinal cord |level |30ug    |SCRaso |  3|  4| -4.1376927| 2.552372| 0.035| 0.0525|ns           |   130.4963|30ug  , SCRaso  |
# ANOVA and corrected posthoc pairwise comparisons
my_aov <- aov(level ~ treatment * tissue, data = my_data)
my_emmeans <- emmeans(my_aov, pairwise ~ treatment | tissue) # emmeans 'adjust' argument is giving odd results, with some p values getting more significant
my_pwc <- as.data.frame(my_emmeans$contrasts) %>%
  dplyr::mutate(group1 = str_split_fixed(contrast, " - ", 2)[, 1],
                group2 = str_split_fixed(contrast, " - ", 2)[, 2],
                p.adj = p.adjust(p.value, method = padj_method),
                p.adj.signif = symnum(p.adj, corr = FALSE, na = FALSE, legend = FALSE,
                                      cutpoints = c(0, 0.001, 0.01, 0.05, 1), 
                                      symbols = c("***", "**", "*", "ns"))) %>%
  relocate(c(group1, group2), .after = "contrast") %>%
  relocate(c(p.adj, p.adj.signif), .after = "p.value")
if (hide_ns) {
  my_pwc_plot <- my_pwc %>%
    dplyr::filter(p.adj < 0.05)
} else {
  my_pwc_plot <- my_pwc
}
y_max <- my_data %>%
  group_by(tissue) %>%
  dplyr::summarise(ymax = max(level, na.rm = TRUE))
if (nrow(my_pwc_plot) > 0) {
  my_pwc_plot <- my_pwc_plot %>%
    dplyr::left_join(y_max, by = join_by(tissue)) %>%
    group_by(tissue) %>%
    dplyr::mutate(y.position = ymax + seq(0.1 * unique(ymax), 0.2 * unique(ymax), length.out = n())) %>%
    ungroup()
} else {
  my_pwc_plot[1,] <- NA
  my_pwc_plot$y.position <- NA
}

# Display ANOVA
print(kable(my_pwc, caption = "ANOVA with corrected posthoc pairwise comparisons"))
## 
## 
## Table: ANOVA with corrected posthoc pairwise comparisons
## 
## |contrast         |group1  |group2 |tissue      |   estimate|       SE| df|    t.ratio|   p.value|     p.adj|p.adj.signif |
## |:----------------|:-------|:------|:-----------|----------:|--------:|--:|----------:|---------:|---------:|:------------|
## |Vehicle - 30ug   |Vehicle |30ug   |Cortex      |  43.011530| 8.026822| 16|  5.3584755| 0.0001792| 0.0002688|***          |
## |Vehicle - SCRaso |Vehicle |SCRaso |Cortex      |  -2.871493| 7.431393| 16| -0.3864003| 0.9213944| 0.9213944|ns           |
## |30ug - SCRaso    |30ug    |SCRaso |Cortex      | -45.883023| 8.026822| 16| -5.7162127| 0.0000896| 0.0001792|***          |
## |Vehicle - 30ug   |Vehicle |30ug   |Spinal cord |  64.887461| 8.026822| 16|  8.0838294| 0.0000014| 0.0000083|***          |
## |Vehicle - SCRaso |Vehicle |SCRaso |Spinal cord |   7.217937| 7.431393| 16|  0.9712764| 0.6047812| 0.7257374|ns           |
## |30ug - SCRaso    |30ug    |SCRaso |Spinal cord | -57.669524| 8.026822| 16| -7.1846021| 0.0000062| 0.0000186|***          |
# Plot
my_plot <-
  ggplot(my_data, aes(x = treatment, y = level)) +
  facet_wrap(vars(tissue), nrow = 1, scales = "free_x") +
  stat_summary(aes(fill = treatment),
               # fun = mean, colour = "black", geom = "bar", alpha = 0.8) +
               fun = mean, colour = "black", geom = "bar") +
  stat_summary(fun.data = mean_se,  
               geom = "errorbar", width = 0.2) +
  geom_jitter(aes(fill = treatment),
              width = 0.2, height = 0, size = point_size, colour = "black", shape = 21) +
  geom_hline(yintercept = 0, colour = "darkgrey") +
  geom_hline(yintercept = 100, colour = "darkgrey", linetype = "dashed") +
  # stat_pvalue_manual(ttest, hide.ns = hide_ns, label = "p.adj.signif")
  stat_pvalue_manual(my_pwc_plot, label = "p.adj.signif", hide.ns = hide_ns, tip.length = 0.01) +
  scale_x_discrete(labels = treatment_ionis_rename,
                   guide = guide_axis(angle = 45)) +
  scale_y_continuous(breaks = pretty_breaks()) +
  scale_colour_manual(values = treatment_ionis_colour) +
  scale_fill_manual(values = treatment_ionis_colour) +
  labs(x = "Treatment",
       y = "Relative expression (%)",
       colour = "Treatment",
       fill = "Treatment") +
  theme_minimal() +
  theme(text = element_text(size = element_text_size),
        legend.position = "none")

# Display plot
print(my_plot)

# Export experiment 2 plot
for (output_format in output_formats) {
  ggsave(filename = file.path(out_dir, "in_vivo_titration", paste0("in_vivo_titration-experiment2_custom.", output_format)),
         plot = my_plot,
         height = 5, width = 7)
}

ASO screen

Fit models

# Experiment ID
experiment_id <- names(MSH3aso_data)[[12]]
experiment_id
## [1] "aso_screen"
# Make results directory
dir.create(file.path(out_dir, experiment_id),
           recursive = TRUE)

# Format data
plot_data <- MSH3aso_data[[experiment_id]] %>%
  pivot_longer(!log_conc_nM,
               names_to = "label",
               values_to = "level") %>%
  separate_wider_delim(label, "_", names = c("treatment", "metric")) %>%
  pivot_wider(names_from = "metric",
              values_from = "level") %>%
  dplyr::mutate(log_conc_nM = as.numeric(log_conc_nM),
                level = as.numeric(level)) %>%
  dplyr::mutate(conc_nM = exp(log_conc_nM)) %>%
  relocate(conc_nM) %>%
  dplyr::filter(!is.na(level))


# Plot using each formula
my_fits <- lapply(model_names, plot_model_function,
                  models = my_models,
                  data = plot_data,
                  xvar = "conc_nM",
                  xlabel = "MSH3 ASO concentration (nM)",
                  yvar = "level",
                  ylabel = "Relative expression (%)",
                  group_var = "treatment",
                  group_label = "Treatment",
                  log_increment = log_increment,
                  plot_se = FALSE)

# Plot model fits
my_plots <- lapply(my_fits, function(my_fit) {
  my_fit$my_plot
})
wrap_plots(my_plots) +
  plot_annotation(title = "Regression models",
                  theme = theme(plot.title = element_text(hjust = 0.5, face = "bold")))

# r2 table
my_r2 <- lapply(my_fits, function(my_fit) {
  my_fit[["r2"]]
})
r2_table <- data.frame(model_name = names(my_r2),
                       r2 = unlist(my_r2)) %>%
  tibble::rownames_to_column("model_number") %>%
  arrange(-r2)

# Display r2 table
kbl(r2_table, caption = "R-squared") %>%
  kable_styling(bootstrap_options = "striped")
R-squared
model_number model_name r2
4 poly3 1.0000000
3 poly2 0.9906616
2 log 0.9801475
1 linear 0.9579879
5 exponential 0.8681231
6 exponential 0.8681231

Plot

# Plot
my_plot <-
  ggplot(plot_data, aes(x = conc_nM, y = level, colour = treatment, fill = treatment)) +
  geom_point(size = point_size) +
  geom_errorbar(aes(ymin = level - sd, 
                    ymax = level + sd),
                width = 0.25) +
  stat_smooth(method = lm, 
              formula = y ~ log(x),
              fullrange = FALSE, alpha = 0.2, se = FALSE) +
  scale_x_continuous(breaks = pretty_breaks()) +
  scale_y_continuous(trans = "log2",
                     breaks = c(0.1, 1, 2, 5, 10, 25, 50, 100)) +
  labs(x = "MSH3 ASO dose (nM)",
       y = "Relative expression (%)",
       colour = "ASO",
       fill = "ASO") +
  theme_bw() +
  theme(text = element_text(size = element_text_size))
my_plot

# Export
for (output_format in output_formats) {
  ggsave(filename = file.path(out_dir, experiment_id, paste0(experiment_id, ".", output_format)),
         plot = my_plot)
}
## Saving 7 x 5 in image
# Plot without legend
my_plot <- my_plot +
  theme(legend.position = "none")
my_plot

# Export
for (output_format in output_formats) {
  ggsave(filename = file.path(out_dir, experiment_id, paste0(experiment_id, "_nolegend.", output_format)),
         plot = my_plot,
         height = 6, width = 6)
}

Viability

# Experiment ID
experiment_id <- names(MSH3aso_data)[[13]]
experiment_id
## [1] "viability"
# Make results directory
dir.create(file.path(out_dir, experiment_id),
           recursive = TRUE)

# Format data
plot_data <- MSH3aso_data[[experiment_id]] %>%
  dplyr::select(-metric) %>%
  pivot_longer(!treatment_duration_weeks,
               names_to = "label",
               values_to = "level") %>%
  separate_wider_delim(label, "_", names = c("treatment", "rep")) %>%
  dplyr::mutate(treatment_duration_weeks = as.numeric(treatment_duration_weeks),
                level = as.numeric(level)) %>%
  dplyr::filter(!is.na(level)) %>%
  dplyr::mutate(treatment_duration_weeks = case_when(treatment_duration_weeks == 9 ~ "9 weeks",
                                                     treatment_duration_weeks == 15 ~ "15 weeks",
                                                     TRUE ~ NA))

# Relevel variables
plot_data <- plot_data %>%
  dplyr::mutate(treatment = fct_relevel(treatment, intersect(treatment_order, unique(plot_data$treatment))),
                treatment_duration_weeks = fct_relevel(treatment_duration_weeks, c("9 weeks", "15 weeks")))

# t-test with correction
ttest <- 
  plot_data %>%
  group_by(treatment_duration_weeks) %>%
  t_test(level ~ treatment) %>%
  adjust_pvalue(method = padj_method) %>%
  add_significance(cutpoints = c(0, 0.001, 0.01, 0.05, 1),
                   symbols = c("***", "**", "*", "ns")) %>%
  add_y_position() %>%
  ungroup()

# Display t-test
kbl(ttest, caption = "T-test") %>%
  kable_styling(bootstrap_options = "striped")
T-test
treatment_duration_weeks .y. group1 group2 n1 n2 statistic df p p.adj p.adj.signif y.position groups
9 weeks level Vehicle MSH3aso-0.022uM 3 3 -2.2792048 2.000000 0.150 0.4860000 ns 185.7200 Vehicle , MSH3aso-0.022uM
9 weeks level Vehicle MSH3aso-0.26uM 3 3 -2.4600632 2.000000 0.133 0.4860000 ns 202.5629 Vehicle , MSH3aso-0.26uM
9 weeks level Vehicle MSH3aso-3uM 3 3 -2.1686631 2.000000 0.162 0.4860000 ns 219.4057 Vehicle , MSH3aso-3uM
9 weeks level Vehicle SCRaso-3uM 3 3 -2.2294431 2.000000 0.156 0.4860000 ns 236.2486 Vehicle , SCRaso-3uM
9 weeks level Vehicle H2O2-500uM 3 3 4.3821316 2.000000 0.048 0.2400000 ns 253.0914 Vehicle , H2O2-500uM
9 weeks level MSH3aso-0.022uM MSH3aso-0.26uM 3 3 -0.4005964 3.897133 0.710 0.9800000 ns 269.9343 MSH3aso-0.022uM, MSH3aso-0.26uM
9 weeks level MSH3aso-0.022uM MSH3aso-3uM 3 3 -0.2154569 3.866679 0.840 0.9800000 ns 286.7771 MSH3aso-0.022uM, MSH3aso-3uM
9 weeks level MSH3aso-0.022uM SCRaso-3uM 3 3 -0.5169365 3.583734 0.635 0.9800000 ns 303.6200 MSH3aso-0.022uM, SCRaso-3uM
9 weeks level MSH3aso-0.022uM H2O2-500uM 3 3 4.7208869 3.999167 0.009 0.1125000 ns 320.4629 MSH3aso-0.022uM, H2O2-500uM
9 weeks level MSH3aso-0.26uM MSH3aso-3uM 3 3 0.1668522 3.997708 0.876 0.9800000 ns 337.3057 MSH3aso-0.26uM, MSH3aso-3uM
9 weeks level MSH3aso-0.26uM SCRaso-3uM 3 3 -0.1521099 3.862324 0.887 0.9800000 ns 354.1486 MSH3aso-0.26uM, SCRaso-3uM
9 weeks level MSH3aso-0.26uM H2O2-500uM 3 3 4.7236434 3.913832 0.010 0.1125000 ns 370.9914 MSH3aso-0.26uM, H2O2-500uM
9 weeks level MSH3aso-3uM SCRaso-3uM 3 3 -0.3012376 3.893224 0.779 0.9800000 ns 387.8343 MSH3aso-3uM, SCRaso-3uM
9 weeks level MSH3aso-3uM H2O2-500uM 3 3 4.4800005 3.885439 0.012 0.1125000 ns 404.6771 MSH3aso-3uM, H2O2-500uM
9 weeks level SCRaso-3uM H2O2-500uM 3 3 4.3568231 3.611482 0.015 0.1125000 ns 421.5200 SCRaso-3uM, H2O2-500uM
15 weeks level Vehicle MSH3aso-0.022uM 2 2 -1.1058058 1.000277 0.468 0.9360000 ns 205.7200 Vehicle , MSH3aso-0.022uM
15 weeks level Vehicle MSH3aso-0.26uM 2 2 -0.8782982 1.000365 0.541 0.9547059 ns 222.5629 Vehicle , MSH3aso-0.26uM
15 weeks level Vehicle MSH3aso-3uM 2 2 -0.6412665 1.000236 0.637 0.9800000 ns 239.4057 Vehicle , MSH3aso-3uM
15 weeks level Vehicle SCRaso-3uM 2 2 -0.9314195 1.000375 0.523 0.9547059 ns 256.2486 Vehicle , SCRaso-3uM
15 weeks level Vehicle H2O2-500uM 2 2 14.6315099 1.055513 0.038 0.2280000 ns 273.0914 Vehicle , H2O2-500uM
15 weeks level MSH3aso-0.022uM MSH3aso-0.26uM 2 2 0.2573233 1.962780 0.821 0.9800000 ns 289.9343 MSH3aso-0.022uM, MSH3aso-0.26uM
15 weeks level MSH3aso-0.022uM MSH3aso-3uM 2 2 0.2794283 1.987604 0.806 0.9800000 ns 306.7771 MSH3aso-0.022uM, MSH3aso-3uM
15 weeks level MSH3aso-0.022uM SCRaso-3uM 2 2 0.2320504 1.955399 0.839 0.9800000 ns 323.6200 MSH3aso-0.022uM, SCRaso-3uM
15 weeks level MSH3aso-0.022uM H2O2-500uM 2 2 2.1475974 1.009965 0.275 0.6530769 ns 340.4629 MSH3aso-0.022uM, H2O2-500uM
15 weeks level MSH3aso-0.26uM MSH3aso-3uM 2 2 0.0508183 1.912147 0.964 0.9800000 ns 357.3057 MSH3aso-0.26uM, MSH3aso-3uM
15 weeks level MSH3aso-0.26uM SCRaso-3uM 2 2 -0.0288608 1.999630 0.980 0.9800000 ns 374.1486 MSH3aso-0.26uM, SCRaso-3uM
15 weeks level MSH3aso-0.26uM H2O2-500uM 2 2 2.0742740 1.013148 0.283 0.6530769 ns 390.9914 MSH3aso-0.26uM, H2O2-500uM
15 weeks level MSH3aso-3uM SCRaso-3uM 2 2 -0.0766325 1.901755 0.946 0.9800000 ns 407.8343 MSH3aso-3uM, SCRaso-3uM
15 weeks level MSH3aso-3uM H2O2-500uM 2 2 1.6052854 1.008506 0.353 0.7564286 ns 424.6771 MSH3aso-3uM, H2O2-500uM
15 weeks level SCRaso-3uM H2O2-500uM 2 2 2.1434571 1.013510 0.275 0.6530769 ns 441.5200 SCRaso-3uM, H2O2-500uM
# ANOVA and corrected posthoc pairwise comparisons
my_aov <- aov(level ~ treatment * treatment_duration_weeks, data = plot_data)
my_emmeans <- emmeans(my_aov, pairwise ~ treatment | treatment_duration_weeks) # emmeans 'adjust' argument is giving odd results, with some p values getting more significant
my_pwc <- as.data.frame(my_emmeans$contrasts) %>%
  dplyr::mutate(group1 = str_split_fixed(contrast, " - ", 2)[, 1],
                group2 = str_split_fixed(contrast, " - ", 2)[, 2],
                p.adj = p.adjust(p.value, method = padj_method),
                p.adj.signif = symnum(p.adj, corr = FALSE, na = FALSE, legend = FALSE,
                                      cutpoints = c(0, 0.001, 0.01, 0.05, 1), 
                                      symbols = c("***", "**", "*", "ns"))) %>%
  relocate(c(group1, group2), .after = "contrast") %>%
  relocate(c(p.adj, p.adj.signif), .after = "p.value")
if (hide_ns) {
  my_pwc_plot <- my_pwc %>%
    dplyr::filter(p.adj < 0.05)
} else {
  my_pwc_plot <- my_pwc
}
y_max <- plot_data %>%
  group_by(treatment_duration_weeks) %>%
  dplyr::summarise(ymax = max(level, na.rm = TRUE))
if (nrow(my_pwc_plot) > 0) {
  my_pwc_plot <- my_pwc_plot %>%
    dplyr::left_join(y_max, by = join_by(treatment_duration_weeks)) %>%
    group_by(treatment_duration_weeks) %>%
    dplyr::mutate(y.position = ymax + seq(0.1 * unique(ymax), 0.5 * unique(ymax), length.out = n())) %>%
    ungroup()
} else {
  my_pwc_plot[1,] <- NA
  my_pwc_plot$y.position <- NA
}

# Display ANOVA posthoc
kbl(my_pwc, caption = "ANOVA with corrected posthoc pairwise comparisons") %>%
  kable_styling(bootstrap_options = "striped")
ANOVA with corrected posthoc pairwise comparisons
contrast group1 group2 treatment_duration_weeks estimate SE df t.ratio p.value p.adj p.adj.signif
Vehicle - (MSH3aso-0.022uM) Vehicle (MSH3aso-0.022uM) 9 weeks -27.000000 26.58425 18 -1.0156390 0.9066705 1.0000000 ns
Vehicle - (MSH3aso-0.26uM) Vehicle (MSH3aso-0.26uM) 9 weeks -34.333333 26.58425 18 -1.2914916 0.7857650 1.0000000 ns
Vehicle - (MSH3aso-3uM) Vehicle (MSH3aso-3uM) 9 weeks -31.000000 26.58425 18 -1.1661041 0.8468018 1.0000000 ns
Vehicle - (SCRaso-3uM) Vehicle (SCRaso-3uM) 9 weeks -37.666667 26.58425 18 -1.4168792 0.7170104 1.0000000 ns
Vehicle - (H2O2-500uM) Vehicle (H2O2-500uM) 9 weeks 52.666667 26.58425 18 1.9811231 0.3896084 1.0000000 ns
(MSH3aso-0.022uM) - (MSH3aso-0.26uM) (MSH3aso-0.022uM) (MSH3aso-0.26uM) 9 weeks -7.333333 26.58425 18 -0.2758526 0.9997426 1.0000000 ns
(MSH3aso-0.022uM) - (MSH3aso-3uM) (MSH3aso-0.022uM) (MSH3aso-3uM) 9 weeks -4.000000 26.58425 18 -0.1504650 0.9999871 1.0000000 ns
(MSH3aso-0.022uM) - (SCRaso-3uM) (MSH3aso-0.022uM) (SCRaso-3uM) 9 weeks -10.666667 26.58425 18 -0.4012401 0.9984202 1.0000000 ns
(MSH3aso-0.022uM) - (H2O2-500uM) (MSH3aso-0.022uM) (H2O2-500uM) 9 weeks 79.666667 26.58425 18 2.9967621 0.0711926 0.5339445 ns
(MSH3aso-0.26uM) - (MSH3aso-3uM) (MSH3aso-0.26uM) (MSH3aso-3uM) 9 weeks 3.333333 26.58425 18 0.1253875 0.9999948 1.0000000 ns
(MSH3aso-0.26uM) - (SCRaso-3uM) (MSH3aso-0.26uM) (SCRaso-3uM) 9 weeks -3.333333 26.58425 18 -0.1253875 0.9999948 1.0000000 ns
(MSH3aso-0.26uM) - (H2O2-500uM) (MSH3aso-0.26uM) (H2O2-500uM) 9 weeks 87.000000 26.58425 18 3.2726147 0.0414274 0.5313155 ns
(MSH3aso-3uM) - (SCRaso-3uM) (MSH3aso-3uM) (SCRaso-3uM) 9 weeks -6.666667 26.58425 18 -0.2507751 0.9998387 1.0000000 ns
(MSH3aso-3uM) - (H2O2-500uM) (MSH3aso-3uM) (H2O2-500uM) 9 weeks 83.666667 26.58425 18 3.1472272 0.0531315 0.5313155 ns
(SCRaso-3uM) - (H2O2-500uM) (SCRaso-3uM) (H2O2-500uM) 9 weeks 90.333333 26.58425 18 3.3980023 0.0321778 0.5313155 ns
Vehicle - (MSH3aso-0.022uM) Vehicle (MSH3aso-0.022uM) 15 weeks -47.000000 32.55892 18 -1.4435368 0.7016682 1.0000000 ns
Vehicle - (MSH3aso-0.26uM) Vehicle (MSH3aso-0.26uM) 15 weeks -32.500000 32.55892 18 -0.9981903 0.9125629 1.0000000 ns
Vehicle - (MSH3aso-3uM) Vehicle (MSH3aso-3uM) 15 weeks -29.500000 32.55892 18 -0.9060497 0.9399425 1.0000000 ns
Vehicle - (SCRaso-3uM) Vehicle (SCRaso-3uM) 15 weeks -34.000000 32.55892 18 -1.0442607 0.8965202 1.0000000 ns
Vehicle - (H2O2-500uM) Vehicle (H2O2-500uM) 15 weeks 44.500000 32.55892 18 1.3667529 0.7452448 1.0000000 ns
(MSH3aso-0.022uM) - (MSH3aso-0.26uM) (MSH3aso-0.022uM) (MSH3aso-0.26uM) 15 weeks 14.500000 32.55892 18 0.4453465 0.9974065 1.0000000 ns
(MSH3aso-0.022uM) - (MSH3aso-3uM) (MSH3aso-0.022uM) (MSH3aso-3uM) 15 weeks 17.500000 32.55892 18 0.5374871 0.9937575 1.0000000 ns
(MSH3aso-0.022uM) - (SCRaso-3uM) (MSH3aso-0.022uM) (SCRaso-3uM) 15 weeks 13.000000 32.55892 18 0.3992761 0.9984568 1.0000000 ns
(MSH3aso-0.022uM) - (H2O2-500uM) (MSH3aso-0.022uM) (H2O2-500uM) 15 weeks 91.500000 32.55892 18 2.8102897 0.1012131 0.6072789 ns
(MSH3aso-0.26uM) - (MSH3aso-3uM) (MSH3aso-0.26uM) (MSH3aso-3uM) 15 weeks 3.000000 32.55892 18 0.0921406 0.9999989 1.0000000 ns
(MSH3aso-0.26uM) - (SCRaso-3uM) (MSH3aso-0.26uM) (SCRaso-3uM) 15 weeks -1.500000 32.55892 18 -0.0460703 1.0000000 1.0000000 ns
(MSH3aso-0.26uM) - (H2O2-500uM) (MSH3aso-0.26uM) (H2O2-500uM) 15 weeks 77.000000 32.55892 18 2.3649432 0.2200472 0.9430596 ns
(MSH3aso-3uM) - (SCRaso-3uM) (MSH3aso-3uM) (SCRaso-3uM) 15 weeks -4.500000 32.55892 18 -0.1382110 0.9999915 1.0000000 ns
(MSH3aso-3uM) - (H2O2-500uM) (MSH3aso-3uM) (H2O2-500uM) 15 weeks 74.000000 32.55892 18 2.2728026 0.2548095 0.9555356 ns
(SCRaso-3uM) - (H2O2-500uM) (SCRaso-3uM) (H2O2-500uM) 15 weeks 78.500000 32.55892 18 2.4110136 0.2040760 0.9430596 ns
# Plot
my_plot <-
  ggplot(plot_data, aes(x = treatment, y = level)) +
  facet_wrap(vars(treatment_duration_weeks)) +
  stat_summary(aes(fill = treatment),
               # fun = mean, colour = "black", geom = "bar", alpha = 0.8) +
               fun = mean, colour = "black", geom = "bar") +
  stat_summary(fun.data = mean_se,  
               geom = "errorbar", width = 0.2) +
  geom_jitter(aes(fill = treatment),
              width = 0.2, height = 0, size = point_size, colour = "black", shape = 21) +
  # geom_violin(trim = TRUE) +
  # geom_boxplot(outlier.shape = NA, width = 0.2) +
  # geom_jitter(aes(colour = treatment, fill = treatment),
  #             width = 0.2, height = 0, size = point_size, alpha = 0.75) +
  # stat_summary(fun = mean, geom = "point", shape = 18, size = point_size + 2, show.legend = FALSE) +
  # stat_summary(fun = mean, geom = "text", show.legend = FALSE, vjust = -1, size = point_size + 2,
  #              aes(label = after_stat(round(y, 1)))) +
  # stat_compare_means(comparisons = pairwise_comparisons,
  #                    label = "p.signif",
  #                    hide.ns = hide_ns) +
  # stat_pvalue_manual(ttest, hide.ns = hide_ns, label = "p.adj") +
  stat_pvalue_manual(my_pwc_plot, hide.ns = hide_ns, label = "p.adj.signif", tip.length = 0.01) +
  geom_hline(yintercept = 0, colour = "darkgrey") +
  geom_hline(yintercept = 100, colour = "darkgrey", linetype = "dashed") +
  scale_colour_manual(values = treatment_colour) +
  scale_fill_manual(values = treatment_colour) +
  scale_x_discrete(labels = \(x) scales::label_parse()(treatment_rename_subscript[x]),
                   guide = guide_axis(angle = 55)) +
  scale_y_continuous(breaks = pretty_breaks()) +
  labs(x = "Treatment",
       y = "Cell viability (%)",
       colour = "Treatment",
       fill = "Treatment") +
  theme_bw() +
  theme(text = element_text(size = element_text_size),
        legend.position = "none")
my_plot

# Export
for (output_format in output_formats) {
  ggsave(filename = file.path(out_dir, experiment_id, paste0(experiment_id, ".", output_format)),
         plot = my_plot,
         height = 6, width = 6)
}

Mihai MSH3 normalised

# Experiment ID
experiment_id <- names(MSH3aso_data)[[14]]
experiment_id
## [1] "mihai_MSH3_normalised"
# Make results directory
dir.create(file.path(out_dir, experiment_id),
           recursive = TRUE)

# Format data
plot_data <- MSH3aso_data[[experiment_id]] %>%
  pivot_longer(everything(),
               names_to = "treatment",
               values_to = "MSH3_counts_normalised") %>%
  dplyr::mutate(MSH3_counts_normalised = as.numeric(MSH3_counts_normalised)) %>%
  dplyr::filter(!is.na(MSH3_counts_normalised))

# Relevel variables
plot_data <- plot_data %>%
  dplyr::mutate(treatment = fct_relevel(treatment, intersect(treatment_order, unique(plot_data$treatment))))

# t-test with correction
ttest <- 
  plot_data %>%
  t_test(MSH3_counts_normalised ~ treatment) %>%
  adjust_pvalue(method = padj_method) %>%
  add_significance(cutpoints = c(0, 0.001, 0.01, 0.05, 1),
                   symbols = c("***", "**", "*", "ns")) %>%
  add_y_position() %>%
  ungroup()

# Display t-test
kbl(ttest, caption = "T-test") %>%
  kable_styling(bootstrap_options = "striped")
T-test
.y. group1 group2 n1 n2 statistic df p p.adj p.adj.signif y.position groups
MSH3_counts_normalised Vehicle SCRaso 9 8 2.015853 14.428636 6.3e-02 0.0630000 ns 722.0400 Vehicle, SCRaso
MSH3_counts_normalised Vehicle MSH3aso_1161149 9 9 22.278108 8.165391 0.0e+00 0.0000001 *** 804.3067 Vehicle , MSH3aso_1161149
MSH3_counts_normalised Vehicle MSH3aso_1161173 9 9 20.646141 8.425659 0.0e+00 0.0000001 *** 886.5733 Vehicle , MSH3aso_1161173
MSH3_counts_normalised Vehicle MSH3aso_1161329 9 9 13.425819 15.976367 0.0e+00 0.0000000 *** 968.8400 Vehicle , MSH3aso_1161329
MSH3_counts_normalised SCRaso MSH3aso_1161149 8 9 16.864602 7.111053 5.0e-07 0.0000011 *** 1051.1067 SCRaso , MSH3aso_1161149
MSH3_counts_normalised SCRaso MSH3aso_1161173 8 9 15.483349 7.286119 8.0e-07 0.0000013 *** 1133.3733 SCRaso , MSH3aso_1161173
MSH3_counts_normalised SCRaso MSH3aso_1161329 8 9 10.428518 14.204306 0.0e+00 0.0000001 *** 1215.6400 SCRaso , MSH3aso_1161329
MSH3_counts_normalised MSH3aso_1161149 MSH3aso_1161173 9 9 -7.666037 13.398977 3.0e-06 0.0000042 *** 1297.9067 MSH3aso_1161149, MSH3aso_1161173
MSH3_counts_normalised MSH3aso_1161149 MSH3aso_1161329 9 9 -3.886817 8.178619 4.0e-03 0.0050000 ** 1380.1733 MSH3aso_1161149, MSH3aso_1161329
MSH3_counts_normalised MSH3aso_1161173 MSH3aso_1161329 9 9 -2.343382 8.459657 4.6e-02 0.0511111 ns 1462.4400 MSH3aso_1161173, MSH3aso_1161329
# ANOVA and corrected posthoc pairwise comparisons
my_aov <- aov(MSH3_counts_normalised ~ treatment, data = plot_data)
my_emmeans <- emmeans(my_aov, pairwise ~ treatment) # emmeans 'adjust' argument is giving odd results, with some p values getting more significant
my_pwc <- as.data.frame(my_emmeans$contrasts) %>%
  dplyr::mutate(group1 = str_split_fixed(contrast, " - ", 2)[, 1],
                group2 = str_split_fixed(contrast, " - ", 2)[, 2],
                p.adj = p.adjust(p.value, method = padj_method),
                p.adj.signif = symnum(p.adj, corr = FALSE, na = FALSE, legend = FALSE,
                                      cutpoints = c(0, 0.001, 0.01, 0.05, 1), 
                                      symbols = c("***", "**", "*", "ns"))) %>%
  relocate(c(group1, group2), .after = "contrast") %>%
  relocate(c(p.adj, p.adj.signif), .after = "p.value")
if (hide_ns) {
  my_pwc_plot <- my_pwc %>%
    dplyr::filter(p.adj < 0.05)
} else {
  my_pwc_plot <- my_pwc
}
y_max <- max(plot_data$MSH3_counts_normalised, na.rm = TRUE)
if (nrow(my_pwc_plot) > 0) {
  my_pwc_plot <- my_pwc_plot %>%
    dplyr::mutate(ymax = y_max) %>%
    dplyr::mutate(y.position = ymax + seq(0.1 * unique(ymax), 0.5 * unique(ymax), length.out = n()))
} else {
  my_pwc_plot[1,] <- NA
  my_pwc_plot$y.position <- NA
}

# Display ANOVA posthoc
kbl(my_pwc, caption = "ANOVA with corrected posthoc pairwise comparisons") %>%
  kable_styling(bootstrap_options = "striped")
ANOVA with corrected posthoc pairwise comparisons
contrast group1 group2 estimate SE df t.ratio p.value p.adj p.adj.signif
Vehicle - SCRaso Vehicle SCRaso 68.97222 25.67849 39 2.685992 0.0744751 0.0930938 ns
Vehicle - MSH3aso_1161149 Vehicle MSH3aso_1161149 504.77778 24.91179 39 20.262605 0.0000000 0.0000000 ***
Vehicle - MSH3aso_1161173 Vehicle MSH3aso_1161173 471.55556 24.91179 39 18.929010 0.0000000 0.0000000 ***
Vehicle - MSH3aso_1161329 Vehicle MSH3aso_1161329 420.00000 24.91179 39 16.859486 0.0000000 0.0000000 ***
SCRaso - MSH3aso_1161149 SCRaso MSH3aso_1161149 435.80556 25.67849 39 16.971622 0.0000000 0.0000000 ***
SCRaso - MSH3aso_1161173 SCRaso MSH3aso_1161173 402.58333 25.67849 39 15.677845 0.0000000 0.0000000 ***
SCRaso - MSH3aso_1161329 SCRaso MSH3aso_1161329 351.02778 25.67849 39 13.670112 0.0000000 0.0000000 ***
MSH3aso_1161149 - MSH3aso_1161173 MSH3aso_1161149 MSH3aso_1161173 -33.22222 24.91179 39 -1.333594 0.6723546 0.6723546 ns
MSH3aso_1161149 - MSH3aso_1161329 MSH3aso_1161149 MSH3aso_1161329 -84.77778 24.91179 39 -3.403119 0.0127230 0.0181758
MSH3aso_1161173 - MSH3aso_1161329 MSH3aso_1161173 MSH3aso_1161329 -51.55556 24.91179 39 -2.069524 0.2536191 0.2817990 ns
# Plot
my_plot <-
  ggplot(plot_data, aes(x = treatment, y = MSH3_counts_normalised)) +
  # stat_summary(fun = mean, geom = "bar", colour = "black", alpha = 0.8) +
  stat_summary(fun = mean, geom = "bar", colour = "black") +
  stat_summary(fun.data = mean_se, geom = "errorbar", width = 0.2) +
  geom_jitter(width = 0.2, height = 0, size = point_size, alpha = 0.6) +
  stat_pvalue_manual(my_pwc_plot, hide.ns = hide_ns, label = "p.adj.signif", tip.length = 0.01) +
  # geom_hline(yintercept = 0, colour = "darkgrey") +
  scale_x_discrete(labels = treatment_rename,
                   guide = guide_axis(angle = 45)) +
  scale_y_continuous(breaks = pretty_breaks()) +
  labs(x = "Treatment",
       y = "Normalised MSH3 counts") +
  theme_bw() +
  theme(text = element_text(size = element_text_size),
        legend.position = "none")
my_plot

# Export
for (output_format in output_formats) {
  ggsave(filename = file.path(out_dir, experiment_id, paste0(experiment_id, ".", output_format)),
         plot = my_plot,
         height = 6, width = 6)
}


# Plot in colour
my_plot <-
  ggplot(plot_data, aes(x = treatment, y = MSH3_counts_normalised)) +
  stat_summary(aes(fill = treatment),
               # fun = mean, geom = "bar", colour = "black", alpha = 0.8) +
               fun = mean, geom = "bar", colour = "black") +
  stat_summary(fun.data = mean_se, geom = "errorbar", width = 0.2) +
  geom_jitter(aes(fill = treatment),
              shape = 21, width = 0.2, height = 0, size = point_size, alpha = 0.8) +
  stat_pvalue_manual(my_pwc_plot, hide.ns = hide_ns, label = "p.adj.signif", tip.length = 0.01) +
  # geom_hline(yintercept = 0, colour = "darkgrey") +
  scale_fill_manual(values = treatment_colour,
                    labels = treatment_rename) +
  scale_x_discrete(labels = treatment_rename,
                   guide = guide_axis(angle = 45)) +
  scale_y_continuous(breaks = pretty_breaks()) +
  labs(x = "Treatment",
       y = "Normalised MSH3 counts") +
  theme_bw() +
  theme(text = element_text(size = element_text_size),
        legend.position = "none")
my_plot

# Export
for (output_format in output_formats) {
  ggsave(filename = file.path(out_dir, experiment_id, paste0(experiment_id, "-colour.", output_format)),
         plot = my_plot,
         height = 6, width = 6)
}

MSH3 protein vs instability

Format data

# Experiment ID
experiment_id <- names(MSH3aso_data)[[15]]
experiment_id
## [1] "MSH3_western_instab_comparison"
# Make results directory
dir.create(file.path(out_dir, experiment_id),
           recursive = TRUE)

# Format data
plot_data <- MSH3aso_data[[experiment_id]] %>%
  dplyr::mutate(clone = "") %>%
  relocate(clone, .after = "genotype")
  
# Add MSH3ko protein level
MSH3ko_protein <- data.frame("experiment_id" = "MSH3ko_CRISPRwt",
                             "differentiation_id" = c("NI1908", "NI1708", "NI1808"), 
                             "genotype" = "MSH3ko", 
                             "clone" = c("Cl27", "Cl37", "Cl26"),
                             "condition_eb" = NA, 
                             "treatment" = "untreated",
                             "dose_uM" = NA, 
                             "treatment_dose" = "untreated", 
                             "group_juliet" = "MSH3ko_untreated",
                             "group_juliet_differentiations" = c("MSH3ko_CRISPRwt-NI1908-Cl27-MSH3ko",
                                                                 "MSH3ko_CRISPRwt-NI1708-Cl37 -MSH3ko",
                                                                 "MSH3ko_CRISPRwt-NI1808-Cl26-MSH3ko"),
                             "MSH3_level_normalised" = 0)
plot_data <- rbind(plot_data, MSH3ko_protein) %>%
  dplyr::mutate(predvar2 = paste0(experiment_id, "-", differentiation_id, "-", clone, "-", genotype))

# Add instability rate
plot_data <- plot_data %>%
  left_join(juliet_slopes %>%
              dplyr::mutate(predvar2 = str_replace_all(predvar2, " ", "")),
            by = join_by(predvar2, group_juliet == predvar1)) %>%
  dplyr::rename(ii_slope = slope)

# Add annotation
plot_data <- plot_data %>%
  left_join(juliet_annotation,
            by = join_by(experiment_id, group_juliet)) %>%
  dplyr::mutate(group_juliet = fct_relevel(group_juliet, unique(juliet_annotation$group_juliet)))

Fit models

# Plot using each formula
my_fits <- lapply(model_names, plot_model_function,
                  models = my_models,
                  data = plot_data,
                  xvar = "MSH3_level_normalised",
                  xlabel = "MSH3 protein level (%)",
                  yvar = "ii_slope",
                  ylabel = "CAG expansion rate (Change in instability index / week)",
                  group_var = NA,
                  group_label = NA,
                  log_increment = log_increment,
                  plot_se = TRUE)

# Plot model fits
my_plots <- lapply(my_fits, function(my_fit) {
  my_fit$my_plot
})
wrap_plots(my_plots) +
  plot_annotation(title = "Regression models",
                  theme = theme(plot.title = element_text(hjust = 0.5, face = "bold")))

# r2 table
my_r2 <- lapply(my_fits, function(my_fit) {
  my_fit[["r2"]]
})
r2_table <- data.frame(model_name = names(my_r2),
                       r2 = unlist(my_r2)) %>%
  tibble::rownames_to_column("model_number") %>%
  arrange(-r2)

# Display r2 table
kbl(r2_table, caption = "R-squared") %>%
  kable_styling(bootstrap_options = "striped")
R-squared
model_number model_name r2
4 poly3 0.8377280
3 poly2 0.8335948
1 linear 0.7820496
2 log 0.5404600
5 exponential 0.4965232
6 exponential 0.4965232

IC50

Dose response model

# Fit dose-response model
my_drm <- drm(ii_slope ~ MSH3_level_normalised, data = plot_data, fct = LL.4())

# Extract IC50 value
ic50 <- ED(my_drm, 50)
## 
## Estimated effective doses
## 
##        Estimate Std. Error
## e:1:50   32.290     16.092
ic50_estimate <- ic50[[1]]
ic50_se <- ic50[[2]]
ic50_upper <- ic50_estimate + 1.96 * ic50_se
ic50_lower <- ic50_estimate - 1.96 * ic50_se

# Find max value for the response variable using the dose response model
max_response_var_drm <- predict(my_drm, newdata = data.frame(MSH3_level_normalised = 100))

# Plot drm model
plot(my_drm, main = "Dose-response curve using drc")
abline(h = 0.5 * max_response_var_drm, col = "blue", lty = 2)
abline(v = ic50_estimate, col = "blue", lty = 2)

# Table of IC50
ic50_table <- data.frame(method = "drc",
                         max_ii_slope = as.numeric(max_response_var_drm),
                         IC50 = ic50_estimate,
                         ci_lower = ic50_lower,
                         ci_upper = ic50_upper)
kbl(ic50_table, caption = "IC50") %>%
  kable_styling(bootstrap_options = "striped")
IC50
method max_ii_slope IC50 ci_lower ci_upper
drc 0.105816 32.29007 0.7495401 63.83061

Interpolation using a dose response model

# Use predicted data from drm curve to calculate the drm curve IC50 by interpolation
# predicted_values_drm <- data.frame(MSH3_level_normalised = seq(0, max(plot_data$MSH3_level_normalised), length.out = 1000))
predicted_values_drm <- data.frame(MSH3_level_normalised = seq(0, max(plot_data$MSH3_level_normalised), length.out = nrow(plot_data)))
predictions <- as.data.frame(predict(my_drm, newdata = predicted_values_drm, interval = "confidence")) %>%
  dplyr::rename(ii_slope = Prediction, 
                ci_lower = Lower, 
                ci_upper = Upper)
predicted_values_drm <- cbind(predicted_values_drm, predictions)
# plot(ii_slope ~ MSH3_level_normalised, data = predicted_values_drm)

# Calculate IC50 using interpolation with the dose response model
ic50_interpolation_drm <- ed_interpolation(data = predicted_values_drm,
                                           response_var = "ii_slope",
                                           pred_var = "MSH3_level_normalised",
                                           formula = as.formula("ii_slope ~ MSH3_level_normalised"),
                                           model_function = "drm",
                                           response_level = 0.5 * max_response_var_drm,
                                           n_bootstrap = 1000,
                                           confidence_level = 0.95)

# Calculate IC100 (response variable 0% of its maximum value)
ic100_interpolation_drm <- ed_interpolation(data = predicted_values_drm,
                                            response_var = "ii_slope",
                                            pred_var = "MSH3_level_normalised",
                                            formula = as.formula("ii_slope ~ MSH3_level_normalised"),
                                            model_function = "drm",
                                            response_level = 0,
                                            n_bootstrap = 1000,
                                            confidence_level = 0.95)

# Plot linear model
plot(ii_slope ~ MSH3_level_normalised, 
     data = ic50_interpolation_drm$predicted_values,
     main = "Dose-response curve using interpolation")
abline(h = 0.5 * max_response_var_drm, col = "blue", lty = 2)
abline(v = ic50_interpolation_drm$ed, col = "blue", lty = 2)
abline(v = ic50_interpolation_drm$ci_lower, col = "lightblue", lty = 2)
abline(v = ic50_interpolation_drm$ci_upper, col = "lightblue", lty = 2)
abline(h = 0, col = "red", lty = 2)
abline(v = ic100_interpolation_drm$ed, col = "red", lty = 2)
abline(v = ic100_interpolation_drm$ci_lower, col = "pink", lty = 2)
abline(v = ic100_interpolation_drm$ci_upper, col = "pink", lty = 2)

# Table of IC50
ed_table <- data.frame(method = "drc_interpolation",
                         max_ii_slope = as.numeric(max_response_var_drm),
                         ed_level = c(50, 100),
                         ED = c(ic50_interpolation_drm$ed, ic100_interpolation_drm$ed),
                         ci_lower = c(ic50_interpolation_drm$ci_lower, ic100_interpolation_drm$ci_lower),
                         ci_upper = c(ic50_interpolation_drm$ci_upper, ic100_interpolation_drm$ci_upper))
kbl(ed_table, caption = "Effective doses (ED)") %>%
  kable_styling(bootstrap_options = "striped")
Effective doses (ED)
method max_ii_slope ed_level ED ci_lower ci_upper
drc_interpolation 0.105816 50 36.81592 36.81562 36.82542
drc_interpolation 0.105816 100 18.13204 18.07901 18.13435

Interpolation using a linear model

# Calculate IC50 by interpolation
my_lm <- lm(ii_slope ~ MSH3_level_normalised, data = plot_data)
max_response_var_lm <- predict(my_lm, newdata = data.frame(MSH3_level_normalised = 100))
ic50_interpolation_lm <- ed_interpolation(data = plot_data,
                                          response_var = "ii_slope",
                                          pred_var = "MSH3_level_normalised",
                                          formula = as.formula("ii_slope ~ MSH3_level_normalised"),
                                          model_function = "lm",
                                          response_level = 0.5 * max_response_var_lm,
                                          n_bootstrap = 1000,
                                          confidence_level = 0.95)

# Calculate IC100 (response variable 0% of its maximum value)
ic100_interpolation_lm <- ed_interpolation(data = plot_data,
                                           response_var = "ii_slope",
                                           pred_var = "MSH3_level_normalised",
                                           formula = as.formula("ii_slope ~ MSH3_level_normalised"),
                                           model_function = "lm",
                                           response_level = 0,
                                           n_bootstrap = 1000,
                                           confidence_level = 0.95)

# Plot linear model
plot(ii_slope ~ MSH3_level_normalised, 
     data = ic50_interpolation_lm$predicted_values,
     main = "Dose-response curve using interpolation")
abline(h = 0.5 * max_response_var_lm, col = "blue", lty = 2)
abline(v = ic50_interpolation_lm$ed, col = "blue", lty = 2)
abline(v = ic50_interpolation_lm$ci_lower, col = "lightblue", lty = 2)
abline(v = ic50_interpolation_lm$ci_upper, col = "lightblue", lty = 2)
abline(h = 0, col = "red", lty = 2)
abline(v = ic100_interpolation_lm$ed, col = "red", lty = 2)
abline(v = ic100_interpolation_lm$ci_lower, col = "pink", lty = 2)
abline(v = ic100_interpolation_lm$ci_upper, col = "pink", lty = 2)

# Table of effective doses
ed_table <- data.frame(method = "lm_interpolation",
                         max_ii_slope = as.numeric(max_response_var_lm),
                         ed_level = c(50, 100),
                         ED = c(ic50_interpolation_lm$ed, ic100_interpolation_lm$ed),
                         ci_lower = c(ic50_interpolation_lm$ci_lower, ic100_interpolation_lm$ci_lower),
                         ci_upper = c(ic50_interpolation_lm$ci_upper, ic100_interpolation_lm$ci_upper))
kbl(ed_table, caption = "Effective doses (ED)") %>%
  kable_styling(bootstrap_options = "striped")
Effective doses (ED)
method max_ii_slope ed_level ED ci_lower ci_upper
lm_interpolation 0.1228149 50 58.71070 39.227615 72.26981
lm_interpolation 0.1228149 100 17.42139 4.648074 27.23633

Dose response curve

# Group colours
my_colours <- juliet_annotation %>%
  dplyr::select(group_juliet, group_colour) %>%
  distinct()
my_colours <- setNames(my_colours$group_colour, my_colours$group_juliet)

# Group rename
my_rename <- juliet_annotation %>%
  dplyr::select(group_juliet, group_rename) %>%
  distinct()
my_rename <- setNames(my_rename$group_rename, my_rename$group_juliet)

# Plot
my_plot <-
  ggplot(plot_data, aes(x = MSH3_level_normalised, y = ii_slope)) +
  geom_point(aes(fill = group_juliet), 
             size = point_size, shape = 21, colour = "black", alpha = 0.8) +
  geom_line(data = ic50_interpolation_drm$predicted_values, 
            aes(x = MSH3_level_normalised, y = ii_slope),
            colour = "black", size = 1) +
  geom_segment(aes(x = -Inf, xend = ic50_interpolation_drm$ed, 
                   y = 0.5 * max_response_var_drm, yend = 0.5 * max_response_var_drm), 
               colour = "blue", linetype = "dashed", size = 0.5, alpha = 0.5) +
  geom_segment(aes(x = ic50_interpolation_drm$ed, xend = ic50_interpolation_drm$ed, 
                   y = -Inf, yend = 0.5 * max_response_var_drm), 
               colour = "blue", linetype = "dashed", size = 0.5, alpha = 0.5) +
  geom_rect(aes(xmin = ic50_interpolation_drm$ci_lower, xmax = ic50_interpolation_drm$ci_upper,
                ymin = -Inf, ymax = 0.5 * max_response_var_drm), 
            fill = "blue", alpha = 0.01, colour = NA) +
  geom_segment(aes(x = -Inf, xend = ic100_interpolation_drm$ed, 
                   y = 0, yend = 0), 
               colour = "red", linetype = "dashed", size = 0.5, alpha = 0.5) +
  geom_segment(aes(x = ic100_interpolation_drm$ed, xend = ic100_interpolation_drm$ed, 
                   y = -Inf, yend = 0), 
               colour = "red", linetype = "dashed", size = 0.5, alpha = 0.5) +
  geom_rect(aes(xmin = ic100_interpolation_drm$ci_lower, xmax = ic100_interpolation_drm$ci_upper,
                ymin = -Inf, ymax = 0), 
            fill = "red", alpha = 0.01, colour = NA) +
  scale_fill_manual(values = my_colours,
                    labels = my_rename) +
  scale_x_continuous(breaks = pretty_breaks()) +
  scale_y_continuous(breaks = pretty_breaks()) +
  labs(x = "MSH3 protein level (%)",
       y = "CAG expansion rate\n(change in instability index / week)",
       fill = "Treatment") +
  theme_bw() +
  theme(text = element_text(size = element_text_size))
my_plot

# Export
for (output_format in output_formats) {
  ggsave(filename = file.path(out_dir, experiment_id, paste0(experiment_id, "_drm.", output_format)),
         plot = my_plot,
         height = 6, width = 9)
}

# Add confidence interval
my_plot <- 
  my_plot +
  geom_ribbon(data = predicted_values_drm,
              aes(ymin = ci_lower, ymax = ci_upper), 
              fill = "darkgrey", alpha = 0.3)
my_plot

# Export
for (output_format in output_formats) {
  ggsave(filename = file.path(out_dir, experiment_id, paste0(experiment_id, "_drm_ci.", output_format)),
         plot = my_plot,
         height = 6, width = 9)
}

Linear model

my_plot <-
  ggplot(plot_data, aes(x = MSH3_level_normalised, y = ii_slope)) +
  geom_point(aes(fill = group_juliet), 
             size = point_size, shape = 21, colour = "black", alpha = 0.8) +
  stat_smooth(method = lm, 
              formula = y ~ x,
              fullrange = TRUE, se = TRUE, colour = "black", fill = "darkgrey", alpha = 0.3) +
  geom_segment(aes(x = -Inf, xend = ic50_interpolation_lm$ed, 
                   y = 0.5 * max_response_var_lm, yend = 0.5 * max_response_var_lm), 
               colour = "blue", linetype = "dashed", size = 0.5, alpha = 0.5) +
  geom_segment(aes(x = ic50_interpolation_lm$ed, xend = ic50_interpolation_lm$ed, 
                   y = -Inf, yend = 0.5 * max_response_var_lm), 
               colour = "blue", linetype = "dashed", size = 0.5, alpha = 0.5) +
  geom_rect(aes(xmin = ic50_interpolation_lm$ci_lower, xmax = ic50_interpolation_lm$ci_upper,
                ymin = -Inf, ymax = 0.5 * max_response_var_lm), 
            fill = "blue", alpha = 0.01, colour = NA) +
  geom_segment(aes(x = -Inf, xend = ic100_interpolation_lm$ed, 
                   y = 0, yend = 0), 
               colour = "red", linetype = "dashed", size = 0.5, alpha = 0.5) +
  geom_segment(aes(x = ic100_interpolation_lm$ed, xend = ic100_interpolation_lm$ed, 
                   y = -Inf, yend = 0), 
               colour = "red", linetype = "dashed", size = 0.5, alpha = 0.5) +
  geom_rect(aes(xmin = ic100_interpolation_lm$ci_lower, xmax = ic100_interpolation_lm$ci_upper,
                ymin = -Inf, ymax = 0), 
            fill = "red", alpha = 0.01, colour = NA) +
  scale_fill_manual(values = my_colours,
                    labels = my_rename) +
  scale_x_continuous(breaks = pretty_breaks()) +
  scale_y_continuous(breaks = pretty_breaks()) +
  labs(x = "MSH3 protein level (%)",
       y = "CAG expansion rate\n(change in instability index / week)",
       fill = "Treatment") +
  theme_bw() +
  theme(text = element_text(size = element_text_size))
my_plot

# Export
for (output_format in output_formats) {
  ggsave(filename = file.path(out_dir, experiment_id, paste0(experiment_id, "_lm.", output_format)),
         plot = my_plot,
         height = 6, width = 9)
}